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SPECTRA  AND  COVARIANCES  FOR  "CLASSICAL"  NONLINEAR  SIGNAL 
PROCESSING  PROBLEMS  INVOLVING  CLASS  A  NON-GAUSSIAN  NOISE 

PART  I.  ANALYTIC  RESULTS  AND  NUMERICAL  EXAMPLES 

1 .  INTRODUCTION 


Non-Gaussian  noise  fields  play  a  critical  role  in  modern 
signal  processing  because  of  the  frequently  dominant  effects  of 
such  noise  and  interference  in  a  wide  variety  of  applications. 
Communication  theory  generally,  and  specifically  telecommunica¬ 
tions,  electromagnetic  and  acoustic  scattering,  man-made  and 
natural  ambient  noise,  optics,  and  underwater  acoustics,  are 
common  areas  of  interest  in  this  respect.  In  the  present  report 
we  are  concerned  primarily  with  underwater  acoustic  noise 
phenomena,  but  the  models  and  results  are  canonical ,  that  is, 
they  take  forms  invariant  to  the  particular  physical  application 
in  question. 

Specif icall> ,  we  are  concerned  with  various  second-order 
statistics  of  non-Gaussian  noise  processes  and  fields  after  they 
have  been  subjected  to  different  types  of  nonlinear  operations, 
such  as  rectification  and  modulation.  A  generic  problem  here  is 
the  passage  of  non-Gaussian  noise  through  a  zero-memory  nonlinear 
(ZMNL)  device.  The  desired  output  statistics  are  typically  the 
mean  (dc),  mean  intensity  (power),  the  covariance  or  correlation 
function,  and  the  associated  spectra.  These  last  include 
wavenumber  spectra  in  the  case  of  noise  fields,  as  well  as  the 
more  general  frequency-wavenumber  spectra  obtained  by  joint 
temporal  and  spatial  Fourier  transformations.  Typical  "class¬ 
ical"  problems  include;  (i)  rectification,  (ii)  determination  of 
output  spectra  and  covariances,  (iii)  calculation  of  ( output ) 
siqnal-to-noise  ratios ,  (iv)  modulation,  (v)  demodulation,  and 
(vi)  special  systems ,  as  for  example,  the  spectriim  analyzer. 

These  and  other  problems  involving  ZMNL  devices  are  described  in 
detail  in  [1;  chapters  5  and  12  -  17].  What  is  new  here  is  the 
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use  of  the  approximate  second-order  probability  density  functions 
and  characteristic  functions  in  the  above  applications  when  the 
noise  processes  are  non-Gaussian. 

A  full  treatment  is  given  in  a  current  study  by  Middleton, 
[2],  which  is  an  expanded  version  of  his  recent  paper  [3],  which 
employs  some  of  the  results  of  the  present  report,  namely,  the 
calculated  covariances  and  spectra.  Here,  we  are  content  to 
summarize  the  pertinent  analytic  results,  the  corresponding 
examples  of  calculated  covariances  and  spectra,  and  the  various 
computational  procedures  associated  with  their  evaluation.  The 
details  of  the  derivations  are  provided  in  [2]  and  [3].  Included 
here,  also,  is  a  selection  of  illustrations  of  the  analytic 
results. 

2.  ANALYTICAL  RESULTS:  A  SUMMARY 

In  the  present  study,  we  address  three  classical  problems 
where  the  goals  are  the  calculation  of  the  covariance  and 
associated  intensity  spectrxim.  Specifically,  we  consider: 

Problem  I.  The  half-wave  v-th  law  rectification  of  Class  A 
noise  fields  and  processes; 

Problem  II.  Phase  modulation  of  a  carrier  by  a  Class  A  noise 
process;  and 

Problem  III.  Frequency  modulation  of  a  carrier  by  a  Class  A 
noise  process. 

Class  A  noise,  as  noted  in  section  3  of  [2],  [3],  is  a 
canonical  form  of  interference  characterized  by  a  coherent 
structure  vis-A-vis  the  (linear)  front-end  stages  of  a  typical 
receiver:  negligible  transients  are  produced  at  the  output  of 
these  stages.  Class  B  noise,  on  the  other  hand,  is  incoherent 
and  highly  impulsive,  such  that  the  front-end  stages  of  the 
receiver  generate  an  output  which  consists  solely  of  (over¬ 
lapping)  transients.  Here,  the  Class  A  models  are  tractable  in 
the  required  second-order  distribution  and  characteristic 
functions,  whereas  the  Class  B  models  are  not  and  must 
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consequently  be  appropriately  approximated  in  second-order;  see 
[4]  and  [5]  for  additional  information.  In  the  present  report, 
we  shall  consider  examples  of  Class  A  noise  only. 

2.1  THE  SECOND-ORDER  CLASS  A  CHARACTERISTIC  FUNCTION 


In  applications  [1]  -  [3],  the  second-order  characteristic 
function,  F2  ( iJ^j ,  iEi2 )  >  plays  a  key  role:  from  it,  we  may  obtain 
the  aforementioned  statistics  of  the  outputs  of  ZMNL  devices, 
spectra  of  angle-modulated  carriers,  and  other  usually  second- 
order  statistics  of  various  nonlinear  operations  arising  in  a 
variety  of  communication  and  measurement  operations. 

(See  [2],  [3]  for  further  discussion.) 

Here,  we  specifically  use  the  approximate  Class  A  noise 
characteristic  function,  F2,  including  an  additive  Gaussian 
component,  given  by 


(2.1) 


where  A  (“A^)  is  the  "overlap"  index,  and  where 


"m+n  ' 


i  o 


A<L^>; 


{2.2b) 


k^/A 


(2.2c) 


and  k. ,  k^  are  the  normalized  covariances  of  the  non-Gauss  and 
Lf  Ca 

Gauss  components,  respectively.  Thus,  qI  ^  1* 
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Here,  p  (=pj.)  is  the  "overlap"  correlation  function 


P(T  '  ) 


1  -  P|t' I  for  P|t'  I  i  1] 

,  p  *=  1/T  , 

0  for  P{t  '  I  >  ij  ® 


(2.3) 


in  which  T  is  the  mean  duration  of  a  typical  noise-signal  of 
®  2  2 

intensity  <B  >/2  =  <L  >.  The  time  delay  t'  is  given  by 


=  T  (=  tj-t^)  , 


(2.3a) 


respectively,  for  space-time  fields  and  received  temporal 

processes.  The  path  delay  (=  |R2  -  Rjl/c^)  accounts  for 

the  time  differential  between  propagation  paths  to  the  points  at 

which  processing  occurs,  cf.  figure  2.1  ff.  Case  A.  The 
•  •  2 

quantities  ^2A  and  Og  are,  respectively,  the  intensity  of  the 
non-Gauss ian  and  Gaussian  components  which  constitute  the  general 
Class  A  model  used  here.  (However,  we  note  that  the  present 
Class  A  model  belongs  either  to  the  strictly  canonical  Class  A 
cases,  where  o'  1  interfering  sources  are  equidistant  from  the 
observer,  or  more  generally,  to  the  much  broader  class  of 
situations  in  which  the  effective  source  distribution  is 
concentrated  in  an  annulus  whose  inner-to-outer  radii  have  a 
ratio  0(1/2)  or  less.  The  former  is  exactly  represented  by  (2.1) 
to  second-order,  while  the  latter  is  approximately  so 
represented,  albeit  a  good  approximation  as  long  as  the 
aforementioned  source  annulus  is  not  too  large.  See 
[5;  section  V,  C],  for  example.  For  an  exact  treatment,  see  also 
[6],  in  an  important  class  of  physical  models.  Finally, 
differentiation  of  F2,  (2.1),  in  the  usual  way,  gives  us  the 
(exact)  covariance  of  the  composite  Class  A  and  Gauss  field, 
namely. 


Ka  +  Kg  , 


(2.4) 


which  in  normalized  form  is 
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^A+G^  ) 


1  +  r' 


(2.4a) 


In  practice,  A  is  usually  less  than  unity,  say  0(0.1  -  0.3) 

mi+m2+n 

typically,  so  that  only  a  comparatively  few  terms  in  A  are 

needed  for  numerical  evaluation  of  (2.1)  and  the  statistical 
quantities  derived  from  it,  cf.  section  2.2  ff.  Note  that  when 
P|t1  i.  1,  p  =  0,  and  r'=  0,  we  get 


m=0  '  ^  n=0 


=  F^(iJ;^)^  '  (2.4b) 

as  expected:  there  is  now  no  correlation  between  process  samples. 
With  a  Gaussian  component,  these  will  be  correlated,  of  course, 
unless  |t|  ^  so  that  kg  0,  cf.  (2.2c). 

2.2  PROBLEM  I:  HALF-WAVE  v-TH  LAW  RECTIFICATION 

(STATIONARY  AND  HOMOGENEOUS  FIELDS) 


Here  we  consider  the  problem  of  obtaining  the  second-order 
(second-moment)  statistics,  M^,  of  a  sampled  noise  field,  a(R,t), 
after  passage  through  a  ZMNL  device,  g,  when  the  noise  is 
generally  non-Gaussian.  Various  processing  configurations  are 
possible.  We  show  two  in  figure  2.1,  below.  Analytically,  we 
have,  for  stationary  and  homogeneous  inputs  [1;  section  2.3-2] 


My(AR,T) 


g(xi)g(x2> 


(2ii) 


I  J  fiUj)  fCUj) 


^1^2 


X  F2(i?;j,iJ;2;AR,T  d^j  di;2  *  Yi  y2  » 


(2.5) 
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where  AR  «  R2-R2/  t  -  and  f(i£;)  is  the  Fourier  transform 

of  the  ZMNL  device  with  ■  y(Rj/tj),  etc.  In  the  present 
cases,  we  have  specifically 

f(i^)  -  p  r(v+l)/{i^)'''^^  ,  V  >  -1  ,  (2.6) 

for  these  half-wave  v-th  law  rectifiers  [1;  (2.101a,b)]. 


SPATIAL  SAMPLING  -*•  •*-  SIGNAL  PROCESSING  - ► 


2 -SENSOR  ARRAY  ZMNL 


Figure  2.1  A.  Two-point  sensor  array  (R2)  giving  sampled  field 
at  two  space-time  points.  B.  A  general  array  (R)  (preformed 
beam),  converting  the  field  a(R,t)  into  a  single  (time)  process 
x(tj).  Both  are  followed  by  ZMNL  devices,  delays,  and  averaging, 
as  indicated  schematically. 


6 


TR  8887 


For  the  Class  A  non-Gaussian  noise  inputs  of  section  2.1 
above,  we  find  that  the  (normalized)  second-moment  for  the 
resulting  rectified  field  is  now 


where  we  have  further  postulated  the  noise  field  to  be  isotropic, 
AR  |ar|,  and  where  specifically. 


J^(Y|mj,m2,n) 


(2.7a) 


Y  s 
a 


fmj^+n 


i  rinr+rT 

-  +  r'  -4 — 


+  r' 


;  a  =  (m^,m2,n)  ,  |Y^|  i  1  . 

(2.7b) 


Specifically,  also,  we  have  the  following  normalized  forms 


My  s  My/2^^2'’/4n  ;  f '  s  pT '  ,  p  -  l/T^  ,  cf.  (2.3)  , 

I  I 

AR  s  AR/Aj^  ,  Aj^  =  correlation  distance,  AR  =  |R2-R2^|  .  (2.8) 


For  numerical  results,  we  select  the  following  models  for  the 
space-time  covariance  functions  of  the  isotropic  and  stationary 
non-Gaussian  and  Gaussian  components  of  the  input  noise  field: 

E  exp(-AR^/A^  -  |(AUj^T'/p)^]  «  exp(-^^  -  |(AUj^f')^]  ;  (2.9a) 

kg  E  expl-AR^/A^  -  i(Au)gT'/P)^)  -  exp(-AR2(Aj^/Ag)2  -  ^{AWgf')^)  i 
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AWg  =  AcOj,/ P  . 


(2.9b) 


Here,  A^  is  a  correlation  distance,  and  Au^^,  AWg  are  angular 
frequency  spreads  associated  with  the  respective  non-Gaussian  and 
Gaussian  components  of  the  input  field.  Note  that  if  we  define 
the  correlation  distance  A^  as  that  where  *  1/e  (f'  *  0), 
then  A^  *  etc. 

For  the  special  cases  of  v  considered  here,  we  also  observe 
(from  [1;  (A. 1-39)])  that  may  be  expressed  in  closed  form: 

Bq(Y)  =  n  +  2  arcsin(Y)  ,  (2.10a) 


Bi(Y) 


arcsin( Y ) 


(>- 


(9.10b) 


B2(Y)  =  (5  +  (f  arcsin(Y)]  +  |y(i  -  Y^j  .  (2.10c) 


2.2-1  GAUSS  PROCESSES  ALONE  (A=0) 


When  only  a  Gauss  noise  field  is  originally  present,  that  is, 
A  =  0,  for  example,  82^  *  0,  (2.7)  reduces  to  the  classical 
result  [1;  page  541,  (13.4a)]: 

(2.11) 

For  comparison  with  the  non-Gaussian  cases  (A>0),  we  choose  to 
have  equal  input  noise  intensities.  This  means  that 

'^A=0  “  ^G  “2A  “  ®2A^^  +  r')  » 

so  that 


/N 

M  ,  o 
y  A=0 


(1  +  r')'’  B^ 


a=0 


(2.12) 


A  A 

and  My  is  then  to  be  compared  with  M^,  A  >  0.  When  r'  is  small, 
as  is  usually  the  case,  we  can  often  replace  (1  +  r')'*  by  unity. 
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At  this  point,  following  figure  2.1,  we  distinguish  two 
classes  of  operation;  (A),  where  a  pair  of  point  sensors  is  used 
to  sample  the  noise  field  and  we  wish  to  consider  both  the  space 
and  temporal  correlations  of  the  sampled  field  at  1  e  two  points 
(Rl,ti),  (R2,t2);  and  {B),  where  the  space-time  field  is 
cci.verted  into  a  random  process,  x(t),  by  the  beamforming  array 
(R),  with  an  associated  directionality  embodied  in  the  resultant 
beam  (vide  [7;  sections  IV  B  and  VI  A]). 


2.2-2  CASE  B,  FIGURE  2.1 

Let  us  consider"  the  simpler  case  (Case  B)  of  the  time  process 
first,  cf.  (B).  For  this,  we  set  AR  =  0  formally  in  (2.7)  et 

A 

seq.  above,  since  x  =  R  a(R,t)  here  and  t'  =  t  = 

(2.3a).  See  also  [3;  (3.2)  et  seq.  and  (3.11a)].  Then  our  ad 

hoc  illustrative  models  of  the  process  covariances  kj^,  kg,  are, 
from  (2.9a,b),  at  once* 

kh  *  *  exp(-  ^(AWj^f/p)^]  «  exp(-  •|(AWj^f)^]  ,  (2.13a) 

kg  *  '  exp(-  ^(AUgf/p)^j  «  exp[-  ■|(AWgf)^]  .  (2.13b) 

Accordingly,  (2.7)  reduces  to 

Case  B:  My(0,f)  s  My(f)g  •  (2.7),  with  =  (2.7b), 

and  (2.13a,b)  and  AR  «  0  therein.  (2.14) 

We  note  that  when  |f|  2  1,  p  =  0,  and  My(0,|f|  i.  1)  reduces  to  a 
simpler  relation  [vis-A-vis  (2.7)],  viz.; 


*  A  physically  derived  model  of  kg  and  kj^  may  be  made  from 
[3;  (3.11a)]  with  L*RL,  R«  (2.9)  etc.,  where  L  is  typically 
given  by  [3;  (3.3)],  for  example. 
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^  -9  A 

My(f)B  =  e 


fm. 


A  A  -1  +  r' 

l_ _ _  nij^  i  m2  I  A 


-4  +  T'  \  B 


where  (2.7b)  becomes 


r'  k. 


.r  _  VJ 

I  —  -7 - r-i - T - — I 

am,  T  m_ 


p  =  0 , 


(2.14a) 


(2.14b; 


in  B 


Special  cases  of  interest  are: 

I.  THE  INTENSITY  E(y^):  f  =  0,  p  =  1,  m^ 
(2.14)  reduce  to 


m2  *  0 ,  and  (2.7), 


2  ^  ^  -A 

y  *=Min)„»M(0,0)=B  *e 

^  norm  y  ^  B  y '  '  '  a=:n 


where  now  y. 


(2.15) 

1,  e.g.,  k, (0)  =  1  etc.,  and  B  is  independent 


of  n,  for  example,  for  Y-  =  1/ 

a 


V  a=n 


2n  for  V  =  0 

n  for  V  =  1  ■  ,  cf.  (2.10)  . 
3ii/2  for  y  =  2, 


(2.16a) 


For  general  v,  Y  ■  1,  we  have  (from  [1;  (A.1-34)]) 


®y  a=n  “  2n^r(v+J5)  ,  v  i  0  . 


(2.16b) 


Thus,  (2.15)  becomes 


Inorm  '  *  2''‘’r(v+>i)  e'*  ^  ^  (2  .  r')'’  . 

n»0  (2.17) 
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The  unnormalized  form  is,  from  (2.8), 


with 


-v-1  . 

y  =  M  (0,0)  =  2^^  r(v+Jj)  Hj'’'(A,r')  (2.18) 

j  1  jjT 


)  =  e  ^ 


n=0 


1  4-  V 


for  V  *  0  , 
for  V  *  1  , 


1/A  +  ( 1+r' )  for  V  =  2  . 


For  other  values  of  v  (>0),  we  must  evaluate  numerically 


(2.18a) 

(2.18b) 

(2.18c) 


II.  THE  MEAN  VALUE,  y;  |f  |  -^  «» 

Now  p  =  0,  n*=0,  Y^=0,  and  (2.7)  reduces  directly,  upon 
use  of  (2.18),  r o 


norm 


v/2 


«y(")B  '  “y'"'”'  *  XI  ^  (S  *  ^'1 

m=0 

.  r='(^]  H<'"'2)(s,r')2  . 

(2.19) 


The  unnormalized  form  of  (2.19)  is,  from  (2.8), 


2'*  2^' 

^  ^2A 

4n 


r^(^)  H['’/2)(A,rM^  /  (2.20) 


and  for  v  even,  we  find,  from  (2.18a,b,c) 


^  1  ,  h}1) 


1  +  r'  ,  ^  i  •  (2.21) 
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7  — ^ 

III.  THE  CONTINUUM  INTENSITY:  y"  -  y 

From  (2.18)  and  (2.20)  we  get  at  once  the  general  result 
for  V  i  0, 


2^  s'* 

“2A 


Tlv+h) 


H 


(V) 


2n 


(2.22) 


which  is  the  generalization  of  (1;  (13.7)],  in  the  classical 

purely  Gaussian  cases,  to  the  present,  dominant  non-Gaussian 

2 

noise  component  22^  (>>  °q) •  In  these  classical  cases,  we  can 
show  at  once  that 


lim  hJ'') 


•2A 


r  ^  (=1 

2-»0  4—r;  ^ 


n=0 


2v 

‘^G 


(=  4''’) 


(2.23) 


2 

where  22^  ®  implies  A  0  and  Bq  -♦  0,  cf.  (2.2b),  so  that 

(2.22)  becomes,  as  expected. 


^  2v  r(V+l5) 

^  ,  2  , 

Gauss  G  1 

4n 

(>  0) 


V  ^  0 


(2.24) 


Figure  13.5  of  [1]  shows  (2.24)  as  a  function  of  rectifier 
law  (V),  as  well  as  (2.18),  (2.20)  in  these  Gaussian  cases.  In 
the  present,  more  general,  situation  of  Class  A  noise,  the 
results  are  more  complex,  as  expected,  with  now  two  additional 
parameters  (A,r'),  descriptive  of  this  much  broader  class  of 
interference. 


12 


TR  8887 


2.2-3:  CASE  A,  FIGURE  2.1 

We  turn  now  to  the  more  general  problem  of  the  covariance  of 
the  Class  A  non-Gauss ian  random  field,  sampled  according  to 
procedure  (A),  shown  schematically  in  figure  2.1  earlier.  Here, 

X  =  a{R,t),  sensed  at  (Rj^,tj^),  (R2,t2),  where  L  *  L,  cf.  (3.3)  in 
[3;  (3.2)].  Equation  (2.7)  applies  here,  with  AR  ^  0  (as  well  as 
for  AR  «  0),  and  we  use  (2.9a,b)  for  our  illustrative  examples, 
which  are  discussed  in  section  3  following.  At  this  point,  we 
recall  from  (2.3a)  that  the  proper  time  delay  to  use  is 
T'  =  T  -  AR/Cq  in  p  =  p(T'),  and  in  some  of  the  structural 
elements  of  the  noise  field  covariances,  cf.  [3;  (3.11b,c)]. 


CASE  I :  f '  =  0 


From  (2.7),  we  have  p  =  1,  =  m2  =  1,  giving 


My(AR,0) 


-A 


B 


a-n 


.•.T 


AR 


,  (2.25) 


n=0 


where  (2.7b)  is  specifically 


Y 


a=n 


n 

A 


+  r' 


kg(AR,0) 


(2.25a) 


For  calculations,  (2.9a,b)  are  used,  with  given  by  (2.7a), 
where  (2.25a)  provides  Y^.  When  AR  «  0,  (2.25)  reduces  to  (2.15) 
et  seq.  for  the  total  intensity  of  the  field  observed  at  Rj^  «  R2. 
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CASE  II:  AR  oo,  |t  ’  I  >  1 

When  AR  -♦  <»,  we  obtain  different  results,  depending  on  t  '  . 
Here  p  »  0,  -+  0,  cf.  (2.9a,b)  in  (2.7b),  and  therefore  n  -  0. 

a 

Accordingly,  (2.7)  becomes 

2 

My(“,|T'|  >  1)  =  My(«>,®)  -  My(0,®)  -=  ,  (2.19)  .  (2.26) 

The  fact  that  AR  -»  ®  ensures  that  0,  a  behavior  similar  to 

that  for  Case  (B)  above,  when  we  consider  the  purely  Gaussian 
noise  process,  section  2.2-i. 

CASE  III:  AR-»®,  0<|t'|<1 

2 

Here,  p  >  0  while  Y^  -»  0,  so  that  B  ,  (2.7b),  becomes  F  (v+ij) 

O  V 

once  more.  The  second  moment  function  (2.7)  is  now 


(2.27) 


which  is  a  minor  simplification  of  (2.7). 


CASE  IV:  AR  -♦  ®,  |t  '  I  *  0 

In  this  special  situation,  where  t  *  AR/c^  ->  ®  in  such  a  way 
that  T '  *  0  and  therefore  p  -  1,  Y^  -  0,  we  obtain  directly  from 

O 

(2.27)  the  comparatively  simple  result. 


My(®,0) 


H 


(V) 


(>  0)  . 


(2.28) 
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2.2-4:  REMARKS 

At  first  glance,  as  AR  -»  <*>,  we  might  expect  M  always  to 

—2  —2  ^ 
reduce  to  y  ,  e.g. ,  -  y  =0  for  the  covariance  of  the 

rectified  space-time  field.  This  is  expectedly  the  case  for  the 

covariance  (and  second-moment)  function  of  the  input  Class  A  and 

Gauss  noise  field  components  a(R,t),  as  we  can  see  directly  from 

{2.9a,b),  or  from  [3;  (3.11b,c)]  for  example,  in  the  physically 

derived  cases.  However,  the  process  or  field  y  =  g(x)  here  is 

the  result  of  a  nonlinear  operation,  cf.  (2.5),  (2.6),  which 

severely  distorts  the  input  waveform  and  generates  all  kinds  of 

modulation  products,  associated  with  the  spatial  as  well  as  the 

temporal  variations  of  the  input  field.  This  accounts  for  the 

_2 

departures  in  Cases  III,  IV  of  M  (®,f')  from  y  ,  while  certainly 
—2  —  ^ 

Mx(<“,T')  ->  X  =0,  (since  x  =  0  initially  here). 

From  the  various  limiting  results  above,  we  see  that 

My(0,0)  >  My(«,0)  and  My(0,0)  >  My(0,oo)  ,  (2.29a) 

and 

My(®,0)  =  My(0,<»)  depending  on  A,  T' ,  and  v  ,  (2.29b) 

with 

M  (0,0)  -  My(0,«)  >  0  ,  cf.  (2.19)  and  (2,20)  , 

_2 

My(”,®)  =  My(0,«»)  =  y  ,  cf.  (2.20)  and  (2.26)  , 
whereas 

_2 

Mx(0,0)  >  M^(0,®)  -  M^(<»,0)  *  X  -  0  .  (2.29e) 

Finally,  we  note  that  (2.11),  (2.12)  apply  here,  also,  for 
the  Gauss-alone  cases,  where  now 

■*  ”1  ®v|a=0  ■*  ®v|a  ■  “  • 


(2.29c) 

(2.29d) 
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2.2-5:  SPECTRA 

The  various  intensity  spectra  associated  with  the  output  of 
the  processor  (cf.  figure  2.1)  are  important  also,  as  they  show 
how  the  energy  in  this  output  is  distributed.  Here,  we  consider 
two  types  of  spectra,  respectively,  for  the  rectified  spatial 
field  (A)  and  for  the  process  (B),  namely  the  wavenumber  and  the 
frequency  spectrum  of  y(R,0)  and  y(0,t).  In  particular, 
wavenumber  spectra  are  useful  in  the  analysis  of  spatially 
distributed  phenomena,  paralleling  the  analysis  of  time-dependent 
phenomena . 


TR  8887 


.  A 

Since  My(<»,0)  is  nonvanishing,  cf.  (2.28),  there  is  a  dc 
component,  or  6-function,  in  the  general  wavenumber  spectrum. 
We  will  use  the  relations 


J  X  J^(kx)  dx  «  &(k  -  0)  ,  k  *  [k^  +  ky]  -=  |k|  ,  k  «  (k,  +  )  , 

0  ^ 

(2.32) 


where  we  must  remember  that  k  is  two-dimensional.  With  v  a 
vector  wavenumber  defined  by 


k  =  2nv  [*  ("0,  +  )]  ,  k  =  2nv  *  2n|v|  , 
and  using  the  relation 

6 (ax  -  b)  =  ^  6 |x  -  “I  for  a  >  0  , 


(2.33a) 


(2.33b) 


we  also  show  that 


6(k  -0)  6(k  -0)  =  6(k-0)  -  - 6(0-0) 

^  2nk  (2n)-^v 


— ^  o(0  -0)  6(0  -0)  ,  0  *=  |0|  ,  V  »  foj  +  oJl 

(2n)‘^  X  y  (  X  y) 


(2.34) 


Applying  (2.32)  -  (2.34),  with 


A 

W2{k,4)y 


A 

+  2n  M 

y 


A 

X  Jq ( kx )  dx 


(2.35a) 


-  «2(''''>ly-cont  ‘(V”’ 

-  W2(<''‘»y-cont  *  "y'"'”’  “V"'  ‘'V'”  ' 


(2.35b) 
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A  A 

which  defines  ^2 { k, 0 the  continuous  portion  of  the 
spectrum  and  shows  the  dc  term  in  k-  or  v-space,  as  convenient. 

A 

It  is  ^2-coTit  which  we  are  concerned  in  the  specific 

numerical  examples  of  section  3  ff. 


II.  FREQUENCY  SPECTRUM 

Here  we  employ  the  Wiener-Khintchine  theorem  [1;  (3.42)]  to 
write  for  the  frequency  spectrum  of  y 

OD  OD 

Wy(f)  “  2  I  My(0,T)  exp(-iwT)  dr  =  J  My(0,f)  cos(wf)  df,(2.36) 

-00  0 

where 

Bq  =  2^/n^j  f  =  Pt;  w  =  2nf;  ca  h  u/p;  f  =  f/p.  (2.36a) 

Accordingly/  we  define  the  normalized  frequency  intensity 
spectrum  of  y  as 


Wy(f)  =  “  J  cos(uf)  df  . 

0 


(2.37) 


cf . 


Again,  there  is  a  dc  component,  since  My(0,<») 
(2.20).  We  have 


(>  0), 


Wy(f)  =  J  [My(0,f  )-My(0,»)]  cos(u)f)  df  +  ^  My(0,«»)  S(f-O),  (2.38) 
0 

since 

00 

J  cos(wx)  dx  “  n  6(<>>-0)  «  ^  &(f-0)  . 

0 
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As  in  the  wavenumber  cases  above  (Case  I),  we  are  concerned  with 
the  continuous  part  of  the  spectrum,  viz. 


(2.39) 


which  is  also  illustrated  numerically  in  section  3  ff. 


III.  WAVENUMBER  FREQUENCY  SPECTRUM 


The  wavenumber  frequency  spectrum  is  defined  by 


W2(k,w)y  =  JJ  My(AR,T)  exp( ik*AR-ia)T  )  d(AR)  dr  ,  (2.40) 

with  b)  «  2nf.  The  associated  wavenumber  spectrum  W2(k,0)  used  in 
(2.30)  is  obtairu;d  from  W2(k,T)|^^Q.  In  normalized  form,  we  have 
for  (2.40),  in  these  isotropic  cases, 

W2(k,w)y  -  (2'’  2^^  A2/(4np))  ^  W2{k,w)y 

CO 

=  JJ  M  (AR,f)  exp( ik-AR-iwf )  d(AR)  df  ; 


00  00 

W2(k/W)y  -  2n  J  J  My(x,f)  JQ(kx)  exp(-iuf)  x  dx  df  . 
0  0 


(2.41) 


The  various  dc  components  are  readily  extracted,  as  in  Cases  I 
and  II  above.  Numerical  examples  of  this  joint  intensity 
spectrum  are  reserved  to  a  possible  subsequent  study.  The 
results  of  section  3  show  the  marginal  spectra  (Cases  I,  II)  of 
this  more  general  situation. 
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2.2-6;  FREQUENCY  AND  PHASE  MODULATION 
BY  CLASS  A  AND  GAUSSIAN  NOISE 

This  is  a  Case  (B)  situation,  cf.  figure  2.1,  where  AR  =  0 
and  we  are  concerned  only  with  the  received  (non-Gaussian)  noise 
process  which  is  used  to  angle-modulate  a  (high  frequency) 
carrier  f^.  For  the  analysis,  see  [3;  section  II]. 

The  general  result  for  the  covariance  of  the  carrier 
modulated  by  Class  A  and  Gauss  noise  is  found  to  be 


^y^^  ^A+G 


^A^  Re[exp(iwQT  -  D^2(T)g  -  A(2-p) 
+  2A(1-p)  expl-D^a^/A])]  , 


(2.42) 


where  now. 


2(r)G 

A 


Cf.  [1;  (4. 2), (14. 14c)], 

k  I 

=  [a^  or  J  <  k|->^)  k(X)  dX 

FM  0 


Dp) 


-  i)  J"(G  or  L)<^>  «  ■  (2.12a) 

0 


Also, 

cf.  [1;  (14.2), 

(14.14c)], 

2(r)G 

-  (“G  ^2A 

)  |k(0)  -  k(T)|,g 

(D^  -  Dp)  {2.42b) 

A 

PM 

and 

o  o 


FM 


J  w^(f)  df/w^  or 
0 


PM 


{2.42c) 
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For  our  numerical  examples,  we  use  the  RC-spectrum  of 
[1;  section  14.1-3],  where  now 

-  exp(-blCl)  ,  ^  '  (2*43) 

and  therefore 

FM:  2(t)^  =  [exp(-blCl)  +  bid  -  1]  ; 

d2  2(T)g  -  [exp(-ld)  +  Id  -  1]  ; 

(^f]a  “  ^2A  •  (2.44a) 


PM;  2(t)^  «  -  exp(-bld)]  ; 

2(T)g  «  H  _  exp(-bld)]  ; 

(^p)a  *  °P  ®2A  '  (2.44b) 


with  b  (>  0)  a  dimensionless  quantity,  as  is  C*  The  quantity 
is  the  bandwidth  of  the  modulating  (Gauss)  noise,  cf. 
(2.43).  Note,  also,  that 


“G  ■’"/‘'“n 


^’(^p)a  -  4  4-  ("pJg  • 


2 

The  quantities  (Pp  p)()  respective  modulation  indexes  for 

FM  and  PM,  cf.  [1;  chapter  14]. 

Finally,  we  have  for  p  in  (2.3),  now  with  AR  ■  0, 


P(T)  -*  p(d 


(2.46) 
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Putting  the  above  (2.43),  (2.44)  in  (2.42),  we  now  specialize 
our  results, 

“  i^o  ^o^^^  cos(fa)QT)  ,  with  k^{0)  »  1  ,  (2.47) 

to  the  normalized  covariance  ^^(t),  respectively,  for  FM  and  PM, 
and  their  associated  spectra.  We  have  for  these  carriers 
modulated  by  a  sum  of  Gaussian  and  Class  A  noise: 

I .  FREQUENCY  MODULATION 


(2.48) 


with  p(T)  given  by  (2.46).  Here,  *  in  (2.42).  Since 

Urn  k  (C)  .  0  , 

there  is  no  dc  in  and  hence  all  the  original  carrier  power 

{~A^/2)  is  distributed  into  the  sideband  continuum  for  this 
highly  nonlinear  modulation,  as  expected  [1;  section  14.1-2]. 

The  associated  intensity  spectrum  for  is  defined  by 

os> 

“<“>a+g|fm  ■  J  dC  ,  a  ■  "SZ  ,  (2.«9) 

0  " 

which  is  determined  by  a  direct  cosine  transform  of  See 

appendix  A. 6  ff. 
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II.  PHASE  MODULATION 


“  ®^p["^'(^p)a  ~  exp(-lcl)]  +2A(1-P)  exp[- 

-  A(2-p)  +  Ap  exp(-  [1  -  exp{-blc|  )  ])  ]  , 

(2.50) 

with  p(C)  again  given  by  (2.46).  We  note  that 

ko(0)pM  =  ^  '  (2.51a) 

as  before;  that  is,  the  total  (normalized)  intensity  is  unity. 
Also 

(2.51b) 

this  is  the  fraction  of  the  power  remaining  in  the  carrier,  so 
that 

•^o^^^PM  ■  ^o^^^PM  '  ^  ■  (2.51b)  ,  (2.51c) 

which  is  the  fraction  of  the  power  distributed  in  the  sideband 
continuiim. 

The  associated  intensity  spectrum  of  the  sideband  continuum 
is  determined  from 

00 

''(“U+G  PM-cont  *  J  [^o(^^PM  "  ^o^"^Pm]  •  (2.5-c) 

0 

See  section  3  ff.  for  examples  and  appendix  A. 5  for  the 
evaluation  methods. 
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Finally,  in  the  equivalent  Gaussian  cases  (Gauss  noise 
modulation  of  equal  intensity  and  basic  spectrum,  e.g., 

^'(4)a  "  -  >■'  4a  <1  ^  '‘g  "  ‘‘l  ' 

we  see  that  (2.48),  (2.50)  reduce  to 

■'o'^'FM-Gauss  ■  ®''p[-(4)a  p'  [«P(-»ltl)  +  b|J.|  -  IJ/b^]  , 

(P^Ia  -  (1  +  r')(4)A  ;  (2-53a) 

‘'o'^'pM-Gauss  ‘  ®’'p[-(4)a  p'  >1  -  exp(-b|Cl)l]  , 

(4)a  -  <1  *  p')(pp)a  - 

with  spectra  obtained  as  before,  from  (2.49)  and  (2.52). 
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3.  NUMERICAL  ILLUSTRATIONS  AND  DISCUSSION 


It  is  convenient  to  discuss  the  general  results,  namely  the 
effects  of  (ZMNL)  nonlinear  rectifiers  on,  and  modulation  by,  a 
mixture  of  Gaussian  and  non-Gaussian  noise  processes  and  fields, 
from  the  specific  numerical  calculations  presented  here  in 
figures  3.1  -  3.10.  These  constitute  a  representative  selection 
from  the  universe  of  possible  parameter  states  [cf.  "Summary  of 
Normalized  Parameters"  and  section  2,  preceding].  This  is  done 
here  on  a  per-figure  basis,  as  noted  below.  In  each  case,  the 
dc  component  is  removed:  only  the  covariance  or  continuous 
spectrum  is  calculated.  We  recall  that  there  are  two  cases  to 
distinguish:  Case  A,  t'=t-AR/Cq,  a  2-element  array;  and  Case  B, 
T'=T=t2-tj^,  a  preformed  beam.  See  figure  2.1  and  (2.3a). 

All  spectra  shown  here  are  normalized  to  have  area  (under  the 
spectrum  level)  of  unity,  i.e.,  the  spectral  normalization  is 
obtained  by  dividing  the  spectrum  by  the  value  of  the  associated 
covariance  at  its  origin.  The  normalization  of  the  covariances 
themselves  is  ootained  by  dividing  by  the  value  at  t^O  or  AR=0. 

I.  GAUSS  NOISE  ALONE 

FIGURE  3.1 

This  figure  shows  the  normalized  temporal  covariance  (AR=0) 
for  both  the  input  and  output  of  a  ZMNL  half-wave  v-th  law  (v20) 
detector,  when  v  *=  0,1,2  and  when  only  Gaussian  noise  (A^O)  is 
applied  to  these  nonlinear  devices.  These  curves  are  based  on 
(2.11)  with  (2.7a),  where  =  kg,  (2.96),  with  Awg  s  AUg/p  *  5 
here.  The  normalization  is  with  respect  to  the  covariance 
maximum;  e.g.,  the  normalized  covariance  shown  in  figure  3.1  is 
obtained  from  [  ( 2 . 11 ) /( 2 . 11  )f  =0  J ,  ^=0.  These  results  apply  for 
both  cases  A,B  of  figure  2.1,  where  now  f'»f,  since  AR=0,  cf. 
(2.3a)  and  remarks . 
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As  expected  (cf.  [1;  chapter  13]),  the  general  nonlinearity 
(2.6),  viO,  contracts  the  covariance,  which  is  equivalent  to 
spreading  the  spectrum  vis-S-vis  the  input  ,  cf.  figure  3.2, 
below.  Moreover,  the  greater  the  distortion  {v=0,2),  usually  the 
greater  are  these  effects.  [See  appendix  A.I.] 

FIGURE  3.2 

This  is  the  same  situation  as  shown  in  figure  3.1,  except 
that  the  normalized  intensity  (frequency)  spectrum  is  calculated 
now  [cf.  section  2.2-5,  Case  II,  (2.39)]  with  My(0,f),  (2.11), 
used  in  (2.39).  Observe  the  greatly  broadened  spectra, 
particularly  at  the  low  spectral  levels,  where  the  greater  spread 
occurs  for  the  "super-clipper",  v=0,  cf.  remarks,  figure  3.1; 
also,  appendix  A. 3. 

FIGURE  3.3 

For  the  same  purely  Gaussian  field  above,  cf.  (2.11)  and 

(2.96),  with  the  spatial  covarirnce  is  calculated,  with 

u 

parameters  *  5^,  using  (2.11)  as  before.  The  normalization 

is  with  respect  to  the  covariance  at  AR*0.  Again,  one  observes 
the  same  kind  of  contraction  in  the  covariance  as  noted  in  figure 
3.1.  [See  appendix  A. 2.] 

FIGURE  3.4 

This  is  the  wavenumber  analogue  of  the  frequency  spectrum  of 
figure  3.2,  now  with  and  is  obtained  from  (2.35a,b)  with 

*  5**.  The  rectification  operation  similarly  spreads  the 

iJ  \9 

wavenumber  spectriun,  with  the  greatest  distortion  (v»0)  yielding 
the  greatest  wavenumber  spread,  as  expected  from  the 
corresponding  contraction  of  the  associated  covariance,  cf. 
figure  3.3  above.  [See  appendix  A. 4.] 
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CLASS  A  PLUS  GAUSS  NOISE 

FIGURE  3.5 

The  temporal  covariance  here  is  given  by  the  general  result 

(2.7),  with  the  associated  relations  (2.7),  (2.8),  (2.9),  wherein 

AR=0,  so  that  T=T'*t2-t^,  as  before,  and  where  g/  (2.7a),  is 

now  given  analytically  by  (2.10)  for  v  «  0,1,2.  Here,  the 

u 

parameter  values  are  AUg  =  AWg/p  *  5^,  as  before,  now  with  A=0.2, 
r^»=10~^,  AwL/p  E  AWj^  *  1  for  the  Class  A  non-Gaussian  noise 
component,  typically. 

Again,  for  the  super-clipper  (V“0),  the  contraction  in  the 
normalized  covariance  is  greatest,  cf.  figure  3.1.  But  the 
contribution  of  the  comparatively  strong  non-Gaussian  component 
exaggerates  this  effect.  [See  appendix  A.l.) 

FIGURE  3.6 

The  corresponding  intensity  (frequency)  spectrum  (AR=0), 
obtained  from  (2.7)  in  (2.39),  however,  shows  a  fine-structure 
not  exhibited  wnen  Gauss  noise  alone  (A=0)  is  applied  to  these 
ZMNL  devices.  The  spectral  levels  for  the  case  v=0,  (A=0)  and 
(A>0),  cf.  figure  3.2  with  figure  3.6,  are  approximately  the 
same,  whereas  the  other  inputs,  cases  v>°l,2,  are  much  elevated  as 
f  becomes  larger,  again  due  to  the  presence  of  the  structured 
Class  A  noise,  when  pT  i  1,  cf.  (2.3);  on  the  average,  the 

5 

original  Class  A  "signals"  are  of  comparatively  short  duration, 
or  spectrally  wide  to  begin  with,  so  that  clipping  further 
spreads  the  spectrum.  [See  appendix  A. 3.] 

FIGURE  3.7 

The  spatial  covariance  when  Class  A  noise  is  added  to  the 
Gaussian  input  shows  analogous  behavior,  cf.  figures  3.3  and  3.5; 
the  covariance  is  compressed  vis-A-vis  the  input,  but  more  so 
than  in  the  Gauss-alone  situations.  Again,  (2.7)  -  (2.10)  are 
employed.  [See  appendix  A. 2.] 
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FIGURE  3.8 

The  corresponding  wavenumber  (intensity)  spectrum  with  Class 
A  noise  and  the  Gaussian  component,  obtained  from  (2.7)  -  (2.10) 
in  (2.39),  is  shown  here.  Comparison  with  figure  3.4  indicates  a 
broader  spectral  input,  due  to  the  non-Gaussian  component,  but  a 
relatively  narrower  output,  although  the  latter  is  still 
noticeably  spread  vis-A-vis  the  original  input.  [See  appendix 
A.4.  ] 

FIGURE  3.9 

Finally,  we  consider  the  angle-modulation  cases  described  in 
section  2.2-6  above,  where  weak  to  strong  angle  modulations 
(p  ^  1  to  50)  by  Class  A  noise,  with  a  weak  (r'=10~^)  Gaussian 
modulation  component,  is  employed. 

For  phase  modulation  by  non-Gaussian  noise,  based  on  (2.50) 
with  (2.44b),  (2.45),  (2.46),  the  resulting  normalized  intensity 
(frequency)  spectra  are  obtained  by  applying  (2.50)  to  (2.52), 
where  f  •=  w/2n;  u  *  (w-Qq)/AWjj,  cf.  (2.49).  Note  the  "spike"  at 
f  ~  0.1,  followed  by  a  variety  of  sidelobes  which  rise  as  the 
phase  modulation  index  pp  increases.  The  spike  is  now  bounded  at 
f  “  0.8,  at  the  -10  dB  level,  when  Pp  »  50.  As  expected,  the 
larger  indexes  (Pp)  produce  broader  spectra.  [See  appendix  A. 5.] 

FIGURE  3.10 

For  frequency  modulation  by  non-Gaussian  noise,  from  (2.48) 
with  (2.49)  and  (2.44a),  the  corresponding  intensity  (frequency) 
spectra  again  exhibit  a  continuous  spike  (f  <  0.1).  With  small 
modulation  indexes  (Pp),  the  spectra  are  less  broad  than  for  the 
larger  indexes,  as  expected.  The  non-Gaussian  noise  component 
dominates  the  spectrum  here.  [See  appendix  A. 6.] 
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EXTENSIONS 

Other  situations  where  the  second-order  Class  A  probability 
density  functions  may  be  applied  are  noted  in  [2]  and  [3].  We 
list  some  of  the  extensions  of  the  analysis  to  the  following 
"classical"  problems: 

1)  The  inclusion  of  representative  signals,  with  Gauss  and 
non-Gauss  (Class  A)  noise,  in  the  problems  already  treated 
here  (section  2); 

2)  The  case  of  the  full-wave  square-law  rectifier,  with  both 
Class  A  and  B  noise,  as  well  as  Gauss  noise; 

3)  The  extension  of  2)  to  include  general  broadband  and 
narrowband  signals; 

4)  The  calculation  of  signal-to-noise  ratios  and  deflection 
criteria,  cf.  [1;  section  5.3-4]. 

5 )  Covariances  and  spectra  for  ZMNL  system  outputs ,  with 
signals  as  well  as  non-Gaassian  noise  inputs; 

b)  The  r61e  of  the  electromagnetic  (or  acoustic)  interference 
(EMI  or  ACl)  scenario,  cf.  [5;  section  2B,5]; 

7)  Evaluation  of  the  large  (FM,PM)  indexes,  or  asymptotically 
Gaussian  cases,  cf.  [12], 

Further  opportunities  to  extend  the  classical  theory  [2], [3], 
now  with  non-Gauss ian  noise  inputs,  are  evident  from  the  examples 
and  methods  described  in  [1;  chapters  5,  12  -  16],  for  instance. 
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FIGURE  3.2  FREQUENCY  (INTENSITY)  SPECTRUM  (FOR  AR-0 ) ;  GAUSS 
NOISE  ONLY;  CF.  (2.39),  USED  WITH  (2.11),  AND  APPENDIX  A. 3 
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FIGURE  3.3  SPATIAL  COVARIANCE  (FOR  f',T»0);  GAUSS 
NOISE  ONLY;  CF.  (2.11),  (2.7b),  AND  APPENDIX  A. 2 
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FIGURE  3.5  TEMPORAL  COVARIANCE  (FOR  AR«=0);  CLASS  A 
AND  GAUSS  NOISE;  CF.  (2.7)-(2.9)  AND  APPENDIX  A.l 
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FIGURE  3.7  SPATIAL  COVARIANCE  (FOR  f',f=0);  CLASS  A 
AND  GAUSS  NOISE;  CF.  (2.7)-(2.10)  WITH  APPENDIX  A. 2 
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FIGURE  3.8  WAVENUMBER  SPECTRUM  (FOR  f',f«=0);  CLASS  A  AND 
GAUSS  NOISE;  CF.  (2.7)-(2.10)  IN  (2.38a,b)  AND  APPENDIX  A. 4 
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NORMALIZED  FREQUENCY 


FIGURE  3.9a  PHASE  MODULATION  (INTENSITY)  SPECTRUM  FOR 
INDEX  ^p=l,2,5,  CLASS  A  AND  GAUSS  NOISE;  CF.  (2.50)  WITH 
(2.44b),  (2.45),  (2.46)  IN  (2.52),  AND  APPENDIX  A. 5 
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FIGURE  3.10a  FREQUENCY  MODULATION  (INTENSITY)  SPECTRUM 
FOR  INDEX  pp=l,2,5,  CLASS  A  AND  GAUSS  NOISE; 

CF.  (2.48)  WITH  (2.49),  (2.44a),  AND  APPENDIX  A. 6 
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NORMALIZED  FREQUENCY 


FIGURE  3.10b  FREQUENCY  MODULATION  (INTENSITY)  SPECTRUM 
FOR  INDEX  /Jp«10,20,50,  CLASS  A  AND  GAUSS  NOISE; 

CF.  (2.48)  WITH  (2.49),  (2.44a),  AND  APPENDIX  A. 6 
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PART  II.  MATHEMATICAL  AND  COMPUTATIONAL  PROCEDURES 

4.  SOME  PROPERTIES  OF  THE  COVARIANCE  FUNCTION 

In  this  section,  we  collect  some  useful  relations  for  the 
covariance  and  auxiliary  functions  encountered  in  the  numerical 
evaluation.  These  are  necessary  for  rapid  computation  of  the 
multiple  series  involved  here  and  also  serve  as  checks  on  the 
numerical  procedures  employed. 

4.1  SIMPLIFICATION  AND  EVALUATION  OF  B^(Y) 

The  function  B^(Y)  is  defined  by  the  following  combination  of 
hypergeometric  functions; 

B^,v,  .r^(^)  r(-  1,  vV 

+  2  +  l]  Y  f(^^-|-^,  I?  Y^)  for  Y^  i  1  .  (4.1) 

For  the  upper  F  function  in  (4.1),  we  have  [1;  (A. 1.39b)] 

V  -  0,  f[o,  0;  i;  Y^)  -  1  ! 

ij 

V  «  1,  f(-  Y^]  *  Y  arcsin(Y)  +  [l  “  Y^]  ; 

V  «  2,  f(-1,  -1;  Y^)  «  1  +  2Y^  ;  (4.2) 

where  arcs in  is  the  principal  value  inverse  sine  function.  On 
the  other  hand,  for  the  latter  F  function  in  (4.1),  we  have 
[1;  (A. 1.39a)  and  (A. 1.39c)] 
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These  three  quantities  can  be  computed  simultaneously  by  the 
following  very  compact  computer  coding  in  BASIC: 

y2=y*Y 

Sq=SQR( 1.-Y2) 

T=ASN(Y) +1.5707963267948966 

B0=T+T 

B1=Y*T+Sq 

B2=( .5+Y2 )*T+1.5*Y*Sq  (4.5) 


Thus,  the  rather  formidable  expression,  above,  for  B^(Y)  can  be 
evaluated  by  the  use  of  just  one  square  root  and  one  arcsin  when 

V  =  0 ,  1 ,  2 . 

The  following  limiting  values,  which  are  obvious,  are  needed 
for  various  special  cases: 
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Bo(0) 

=  n  , 

Bod)  = 

2n  , 

B,  (0) 

=  1  , 

B3(1)  = 

n  / 

B2(0) 

=  n/4  , 

B2(1)  = 

3n/2  , 

B3(0) 

=  1  , 

B3(1)  = 

15n/4  , 

B4(0) 

=  9n/16  , 

B4(1)  = 

105n/8  . 

These  are  special  cases  of 

B^(0)  =  r^j 

fv  +  1] 
12]' 

B^(l)  =  2  r(v  +  i)  ,  (4.8) 

the  latter  following  from  (10;  (15.1.20)]. 

4.2  LIMITING  VALUES  OF  THE  COVARIANCE  FUNCTION 


The  covariance  function  at  normalized  separation  AR  and 
delay  f  is  given  by  (2.7)  as 


W  00 


My(^R,f  ) 


=  exp[-A(2-p) ] 


(Ad  -  P)] 


inj=0  ni2=0 


in  ^  •  ^2  * 


n=0 


n  +  nil 

V/  ^ 

n  +  m. 

i  _L  Tn# 

Z  .  pf 

A  ^  ^A 

A  A 

v/2 


B^(Y)  , 


(4.9) 


where 


p ( f )  *  max { 0 ,  1  -  I f I }  , 


(4.10) 


Y  -  Y(mj^,m2,n) 


V  +  T' 
Ah  AG 


n  +  ra^ 

- s 

T" 

n  +  m«  1 

*  +  r; 

+  r' 

A 

Aj 

A  AJ 

(4.11) 
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,  =  exp  -  i  ^  ^  , 


(4.12) 


^  f  f^Tl2  ,fAw_12  -1 

=  k^(AR,f )  =  exp  -  1-;^  i  f^l  . 


(4.13) 


The  functions  p,  kj^,  kg  can  be  replaced  by  other  functional 
dependencies,  if  desired.  The  function  B^(y)  has  been 
considered  earlier  and  considerably  simplified  for  v  =  0,  1,  2. 

4.3  VALUE  AT  INFINITY 

As  AR  or  f  -»  ±®,  then 

p-*0,  kj^^O,  kg->0,  Y-»0.  (4.14) 

(If  |f I  remains  less  than  1  as  ^  tends  to  infinity,  then  p  does 
not  approach  zero;  this  nuance  has  been  discussed  elsewhere  in 
this  report.)  Then,  it  follows  that 


because  the  sum  on  n  can  be  terminated  with  the  n  =  0  term. 

The  sum  on  m  can  be  effected  in  closed  form,  for  v  =  0,  2,  4, 
etc.,  by  using  the  following  results; 


=  exp(A)  , 


(4.16) 
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00 


m=0 


00 


m=l 


=  A  exp(A)  , 


00  00 


in=0  m=l 


(m  -  1  +  1 ) 


(4.17) 


w 

-r 


Z—  (m  -  2)  ! 

in=2  ni=l 


w 

(m  -  1)!  ~ 


(A  +  A)  exp(A)  .  (4.18) 


There  follows 


for  V  =  O' 

for  V  =  2  •  .  (4.19) 

for  V  =  4 


The  case  for  v  =  1  requires  a  numerical  summation,  once  A  and 

are  specified.  When  these  limiting  values  are  subtracted  from 
the  correlation  function,  we  obtain  the  covariance  function. 
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4.4  VALUE  AT  THE  ORIGIN 

For  AR  =  0,  f  =  0,  then 

p  =  1,  =  1/  ^  '  (4.20) 


and 


My(0,0) 


00 

=  exp(-A)  ^  It  (i  +  ^k] 


B^(l) 


(4.21) 


n=0 


because  the  sums  on  m^^  and  m2  can  be  terminated  with  the  zero 
terms,  thereby  also  leading  to  Y  =  1. 

The  sum  on  n  can  be  accomplished  in  closed  form,  for 
V  =  0,  1,  2,  etc.,  by  using  results  given  earlier.  There  follows 


r2n 


for  V  =  0 


My(0,0) 


+ 


for  V  =  1?  . 
for  V  =  2 


(4.22) 
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PART  III.  APPENDICES  AND  PROGRAMS 

APPENDIX  A.l  -  EVALUATION  OF  COVARIANCE  FUNCTION 
FOR  ZERO  SEPARATION  =  0) 

A  program  for  the  numerical  evaluation  of  covariance 
M  (Z\R,t)  for  AR  =  0  is  contained  in  this  appendix.  Inputs 

y  2  2 

required  of  the  user  are  A,  (AWj^/p)  ,  (AUq/^)  ,  S(f),  N(f), 

in  lines  20  -  70.  Since  we  are  generally  interested  in  values  of 

A  less  than  1,  the  series  for  fty  in  (4.9)  will  not  have  to  be 

taken  to  very  large  values  of  m, ,  m,,  n;  accordingly,  the  values 
k  1  z 

of  {A  /k! 1  are  tabulated  once  in  lines  260  -  300  with  a  tolerance 
of  lE-10  set  in  line  80. 

The  values  of  the  covariance  at  infinity,  as  given  by  (4.19), 

are  computed  and  subtracted  in  lines  220  -  240  and  400  -  420; 

this  is  in  anticipation  of  taking  a  Fourier  transform  of  a 

covariance  function  which  decays  to  zero  for  large  arguments  AR. 

The  functions  B^(Y)  and  My(AR,f)  are  available  in  the  two 

subroutines  starting  at  lines  1010  and  1120,  respectively.  The 

latter  subroutine  actually  calculates  the  covariance  at  general 

nonzero  values  of  both  AR  and  f,  although  we  only  employ  it  for 

AR  «  0  in  this  appendix;  see  lines  10  and  380.  Also,  for  AR  «  0, 

2 

the  parameter  Lg2  ■  (Aj^/Ag)  xs  not  relevant  and,  hence,  is 
entered  as  zero  in  line  380. 

The  exponential  Gaussian  forms  for  k^  and  kg  are  used  in 
lines  1200  and  1210,  while  the  triangular  form  for  p  is  entered 
in  line  1240.  Any  of  these  can  be  replaced,  if  desired,  by  forms 
more  appropriate  to  the  user. 

The  program  is  written  in  BASIC  for  the  Hewlett  Packard  9000 
Computer  Model  520.  The  designation  DOUBLE  denotes  integer 
variables,  not  double  precision.  The  output  from  the  program  is 
stored  in  data  files  AOTO,  AOTl,  A0T2,  for  v  «  0,  1,  2, 
respectively. 
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10 
20 
30 
40 
50 
60 
70 
80 
90 
100 
lie 
120 
130 
140 
150 
160 
170 
1  80 
190 
200 
210 
220 
230 
240 
250 
260 
270 
280 
290 
300 
310 
320 
330 
340 
350 
360 
370 
380 
390 
400 
410 
420 
430 
440 
450 
460 
470 
480 
490 
500 
510 
520 
530 
540 
550 


Rc=0. 
fl=.  2 
Gp= . 00 1 

in  b2=l  . 
Wgb2=25. 
Dtc=. 01 
Htc=200 


DelR- 

fl<sub8) 

GfiMMfl" (subft) 
<DelW<subL)/Beta)‘'2 
<DeUK£ubG)  /Bet  a>-'2 
INCREMENT  IN  Tau^ 
NUMBER  OF  Tau^  VALUES 


Toler ance=l.E-10 

COM  fif  <  0 : 40  > , C  <  0 ; 80  > , Sq  <  0 : 80  > 

COM  DOUBLE  J  !  INTEGER 

LIN  KagC200) ,  Tc  <0;  20O:» ,  FOCO:  200) ,  FI  <0:  200> ,  F2':0:  200) 

DOUBLE  Nic,K  !  INTEGERS 

FOR  K=0  TO  NIC 


Tc=K*Dtc  !  Tau'' 

Rho=MflX(0. , 1 . -RESCTc ) )  !  Rho 

T2=. 5*Tc*Tc 


K1  =EXP<:-W1  b2*T2) 


^'g  =  EVP<-Ugb2*T2) 

Kag<K)  =  <:Rho*K1 +Gp*Kg)/<  1  . +Gp)  !  INPUT  COVARIANCE 

NEXT  K 

A1  =  1..'A  !  R>0  REQUIRED 

F0i  t-if  =  PI 

FI  i  rif  =  FNFl  i  nf<  A,  Gp) 

F2inf=.25*PI*<l.+Gp)*<l.+Gp) 

Af <0)*1 . 

FOR  K=1  TO  40 
J  =  K 

Af  <K)=T  =  Rf  nC-l  )*A/'K  !  A''K/K! 

IF  T<To1erance  THEN  320 
NEXT  K 

PRINT  "40  TERMS  IN  Af<*)" 

FOR  K=0  TO  J*2 
C<K)=T=K*A1+Gp 
Sq<K)=l . /SQR<T) 

NEXT  K 

FOR  K=0  TO  Ntc 

Tc (K)=Tc=K*Dtc  !  Tau- 

CALL  Nyc<Rc,Tc,A,Gp,W1b2,Wgb2,0. , F0<K) , F 1 (K ) , F2 (K ) ) 
NEXT  K 

MAT  F0=F0-<F0inf ) 

MAT  F1=F1-<F1 inf ) 

MAT  F2=F2-<F2i nf ) 

MAT  F0  =  F0/<F0<0)  ) 

MAT  F1=F1/<F1 <0) > 

MAT  F2=F2/<F2<0) ) 

PRINT  "INFINITY; " ; F0 i nf ; F 1 i nf ; F2 i nf 

PRINT  "MININA:  " ; M I N  <  F0  <  * ) ) ; M I N  <  F 1 ( *  > ) ; M I N ( F2  <  * ) ) 

PRINT  "AT  Ntc;  " ; F0 C Nt c ) ; F 1 < Nt c ) ; F2 < Nt c ) 

CREATE  DATA  "AIT1",8 
ASSIGN  #1  TO  "AITl" 

PRINT  i»l|Kag(*) 

CREATE  DATA  "AOT0",8 
ASSIGN  «1  TO  "AOT0" 

PRINT  #ljF0<*) 

CREATE  DATA  "R0T1",8 
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560  hSSIGH  #1  TO  “FlOTl" 

570  PRINT 

580  CREATE  DATA  "A0T2",8 

590  ASSIGN  TO  "A0T2" 

600  PRINT  #1;F2(*) 

610  ASSIGN  #1  TO  * 

620  Tc  ru5ix  =  Dt  c  »Nt  c 

630  GINIT  200/260 

640  PLOTTER  IS  505,"HPGL" 

650  PRINTER  IS  505 

660  LIMIT  PLOTTER  505,0,200,0,260 

670  VIEWPORT  22,85,19,122 

680  WINDOW  0. . 1 . . 0. , 1 . 

690  PRINT  "VSS" 

700  GRID  . 25, . 25 

710  PRINT  "VSSS" 

720  PLOT  Tc<*),Kag<*) 

730  PENUP 

740  PLOT  Tc(*),F0<*) 

750  PENUP 

760  PLOT  Tc<*),Fl<*) 

770  PENUP 

780  PLOT  Tc  (♦) ,  F2<*.'‘ 

790  PENUP 

800  PAUSE 

810  PRINTER  IS  CRT 

820  PLOTTER  505  IS  TERMINATED 

830  END 

840  ! 

850  DEF  FNF 1 1 nf < A , Gp )  !  for  v<=nu)  -  1 

860  Tol=l.E-l8 

870  flg=R*Gp 

880  T=l. 

890  S=SQR<:  1 . +Rg) 

900  FOR  M=2  TO  100 

910  T=T*R/M 

920  P=T*SQR<M+Ag) 

930  S=S+P 

940  IF  P<S*Tol  THEN  970 

950  NEXT  M 

960  PRINT  "lOO  TERMS  IN  FNFlinf" 

970  T=Gp+A*S*S+2 • *SQR< flg ) ♦S 

980  RETURN  EXP<-2.*A)*T 

990  FNEND 

1000  ! 

1010  SUB  Bnu<Y, B0,  B1  ,  B2)  !  Bo<Y>  for  o=0,l,2 

1020  IF  Y>1.  THEN  PRINT  "Y  =  I  +";Y-1. 

1030  IF  Y>1.  THEN  Y=l. 

1040  Y2  =  Y*Y 

1050  Sq=SQR< 1 . -Y2) 

1060  T  =  ASNi:Y>  +  1 . 570796326794  8966 

1070  B0  =  T  +  T 

1080  B1=Y*T+Sq 

1090  B2=< .5+Y2)*T+l .5*Y*Sq 

1100  SUBEND 

1110  ! 


51 


TR  8887 


1120  SUE  Mmc  <Rc  ,  Tc  ,  fl,  Gp,  W1  b2,  llgb2,  Lg2,  S0,  SI ,  S2> 

1130  COM  Fif  C  *  )  ,  C  C  *  ) ,  SqC  *  ) 

1140  COM  DOUBLE  J  !  INTEGER 

1150  RLLOCFlTE  flp  <  0 :  J  )  ,  flp  1  <  0 :  J  > 

1160  DOUBLE  K,M1,M2,U,K1,K2  !  INTEGERS 

1170  h1=1./R  !  8>G  REQUIRED 

1180  T2=.5*Tc*Tc 

1190  R2=Rc*Rc 

1200  K1 =EXP(-R2-W1 b2*T2) 

1210  Kg  =  EXP(:-Lg2*R2-Wgb2*T2) 

1220  fll;=Rl*Kl 

1230  Gk=Gp*Kg 

1240  Rho  =  Ml=lX<0.  ,  1 . -RBSCTc  )  )  !  Rho 

1250  Rhol=l.-Rho 

1260  fip(0)=flpl <0)=Pk=Pk  1  =  1 . 

1270  FOR  K=1  TO  J 

1280  Pk=Pk*Rho 

1290  Pk l=Pk l*Rbo1 

1300  T=fif<K) 

1310  flp<K)=T*Pk 

1320  flpl(K)=T*Pkl 

1330  NEXT  K 

1340  S0tril=Slml=S2tnl=0. 

1350  FOR  f11=0  TO  J 

1360  S0tri2  =  Slffi2  =  S2m2  =  0. 

1370  FOR  M2=0  TO  J 

1380  S0n=Sln=S2n=0. 

1390  FOR  N=0  TO  J 

1400  K1=H+M1 

1410  K2  =  H  +  ri2 

1420  T=flp<N) 

1430  P=C(K1 )*C<K2) 

1440  Y=<N*ftk+Gk )*Sq<Kl )*Sq<K2) 

1450  CALL  Bnu<Y, B0, B1 , B2) 

1460  S0n=S0n+T*B0 

1470  Sln  =  Sln  +  T*SQR<:P)*Bl 

1480  S2n=S2n+T*P*B2 

1490  NEXT  N 

1500  T2=flpl<M2) 

1510  S0tn2  =  S0rri2  +  T2*S0n 

1520  Slr(i2  =  Slm2  +  T2*Sln 

1530  S2iti2  =  S2iri2  +  T2*S2n 

1540  NEXT  M2 

1550  Tl=Flpl<Mi> 

1560  S0ml=S0ff.l+Tl*S0m2 

1570  Slrnl=Slrnl+Tl*Sltri2 

1580  S2ml=S2ml+Tl#S2m2 

1590  NEXT  Ml 

1600  T=EXP<-fl*<2. -Rho) ) 

1610  S0  =  T*S0iiil 

1620  Sl=T*Sltnl 

1630  S2  =  T*S2tiil 

1640  SUBEND 
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APPENDIX  A. 2  -  EVALUATION  OF  COVARIANCE  FUNCTION 
FOR  ZERO  DELAY  (f,f'  =  0) 


A  program  for  the  numerical  evaluation  of  covariance 
M  (AR,t)  for  T,T'  *0  IS  contained  in  this  appendix.  Inputs 

y  ,  2  ^ 

required  of  the  user  are  A,  (Aj^/Ag)  ,  6(AR),  N(AR),  in  lines 

20  -  60.  The  tolerance  for  terminating  the  triple  infinite  sums 
is  set  at  lE-15  in  line  70.  The  output  from  the  program  is 
stored  in  data  files  AORO,  AORl,  AOR2,  for  v  =  0,  1,  2, 
respectively.  Other  relevant  comments  are  made  in  appendix  A.l. 

A  Ak 

The  limit  of  My  at  AR  =  “  (when  f  =0)  is  given  by  the  closed 
form  results 

fn  for  V  =  01 


My(«,0)  =  1  + 


for  V  *  1 >  . 


(A. 2-1) 


+  fi  + 

[H  Ir.  AJ 


for  V  =  2 


These  values  have  been  subtracted  from  M  so  that  we  can  Bessel 

y  Ak 

transform  a  function  which  tends  to  zero  as  AR  -♦  ®. 


Tc=0.  ! 

8=.  2  ! 

Gp=.001  ! 

Lg2=5.  ! 

Drc=.005  ! 

Nrc=900  ! 

Toler an ce=l.E-15 
COM  flf <0; 40) , C<0: 80) , Sq<0: 80) 
COM  DOUBLE  J  ! 


T  au'' 
fi<  subfl ) 

GfiMMA' <subR) 
<DelL/De1G)''2 
INCREMENT  IN  DelR- 
NUMBER  OF  DelR-  VALUES 


INTEGER 


DIN  Rc  <0;  900)  ,  Kag':0;  900)  ,  F0<0;  900)  ,  FI  <0:  900)  ,  F2C0;  900) 


DOUBLE  Nrc,K 
Al=l . /A 
F0i nf=PI 
F 1 i nf = 1 . +Gp 

F2  i  nt'  = .  25*P  I  *  <  <  1  .  +Gp  )  ■*■  <  1 . 
Rf <0)=1 . 

FOR  K.=  l  TO  40 
J  =  K 

flf <K)=T=Af <K-1 )*R/K 
IF  T<Tolerance  THEN  230 
NEXT  K 

PRINT  "40  TERNS  IN  Rf<*)" 
FOR  K=0  TO  J*2 
C<K)=T=K*A1+Gp 
Sq<K)=l . rSQRCT) 

NEXT  K 


+Gp)+Al ) 


INTEGERS 
A>0  REQUIRED 

!  LIMITS  FOR 
!  Rc=i nf i ni ty 
Al)  !  AND  Tc=0 
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270  FOR  K=0  TO  Hrc 

280  Rc  (K)=Rc=K*Drc  t  DeIR'' 

290  R2=Rc*Rc 

3W0  K1=EXP(-R2> 

310  Kg=EXP<-Lg2*R2> 

320  Rho  =  Mfi>i<:0.  ,  1  . -flBS<Tc  )  >  !  Rho 

330  Kag<lO  =  i;Rho*ia +Gp*Kg)/<:  1  . +Gp)  !  INPUT  COVflRIflHCE 

340  CALL  Myc(:Rc,Tc,fl,Gp,0.  ,0.  ,Lg2,F0<;K>,Fl<lO,F2<K)> 

350  NEXT  K 

360  MAT  F0=F0-<F0inf ) 

370  MAT  F1=F1-<F1 inf ) 

380  HAT  F2  =  F2-<:F2inf ) 

390  MAT  F0=F0/<F0(0) ) 

400  NAT  F1=F1/<F1<0)) 

410  MAT  F2  =  F2.'<F2<0)  ) 

420  PRINT  "INFINITY: ‘';F0inf;Flinf;F2inf 

430  PRINT  "MINIMA;  " ; H I N < F0 < ♦ ) ) ; M I H< F 1 < * )  )  ; M I N < F2 ( * >  ) 

44  0  PRINT  "AT  Nrc  ;  "  ;  F0  Nr  c  >  ;  F  1  (  Nr  c  )  ;  F2  (  Nrc  > 

450  CREATE  DATA  "AIR1",33 

460  ASSIGN  #1  TO  "AIRl" 

470  PRINT  iH;Kag<*) 

480  CREATE  DATA  "ROR0",33 

490  ASSIGN  #1  TO  "AOR0" 

500  PRINT  #1;F0<:*) 

510  CREATE  DATA  "A0R1",33 

520  ASSIGN  #1  TO  "AORl" 

530  PRINT  #1|F1<»> 

540  CREATE  DATA  “A0R2",33 

550  ASSIGN  #1  TO  "A0R2" 

560  PRINT  #i;F2<*) 

570  ASSIGN  TO  * 

580  Rc  fnax  =  Drc  *Nrc 

590  GINIT  200/260 

600  PLOTTER  IS  505,"HPGL" 

610  PRINTER  IS  505 

620  LIMIT  PLOTTER  505,0,200,0,260 

630  VIEWPORT  22,85,19,122 

640  WINDOW  0. , 3. , 0. , 1 . 

650  PRINT  "VS5" 

660  GRID  .5, .25 

670  PRINT  "VS36" 

680  PLOT  Rc <*) , KagC*) 

690  PENUP 

700  PLOT  Rc<*),F0(*) 

710  PENUP 

720  PLOT  Rc<*),Fl<*) 

730  PENUP 

740  PLOT  Rc<*),F2<*) 

750  PENUP 

760  PAUSE 

770  PRINTER  IS  CRT 

780  PLOTTER  505  IS  TERMINATED 

790  END 

800  I 

810  SUB  Bnu<Y,B0,Bl,B2)  !  Bu<Y)  for  o=0,l,2 

820  !  SEE  APPENDIX  A. 1 

900  SUEEND 

910  ! 

920  SUB  Myc(Rc,Tc,R,Gp,W1b2,Wgb2,Lg2,S0,Sl,S2> 

930  !  SEE  APPENDIX  A. 1 

1440  SUBEND 
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APPENDIX  A. 3  -  EVALUATION  OF  TEMPORAL  INTENSITY 

SPECTRUM  FOR  ZERO  SEPARATION  ( AR  =  0 ) 


A  program  for  the  numerical  evaluation  of  the  Fourier 
transform  of  covariance  My(0,f)  -  My(“)  is  contained  in  this 
appendix.  Inputs  required  of  the  user  are  listed  in  lines 
10  -  30.  The  data  input,  AOTO  or  AOTl  or  AOT2,  as  generated  by 
means  of  the  program  in  appendix  A.l,  is  injected  by  means  of 
lines  410,  600,  and  790. 

In  order  to  keep  the  FFT  (fast  Fourier  transform)  size,  N  in 
lines  30  and  320,  at  reasonable  values,  the  data  sequence  is 
collapsed,  without  any  loss  of  accuracy,  according  to  the  method 
given  in  [8;  pages  7-8]  and  [9;  pages  13  -  16].  The 
integration  rule  documented  here  is  the  trapezoidal  rule;  this 
procedure  is  very  accurate  and  efficient  and  is  recommended  for 
numerical  Fourier  transforms. 


10 

Nt.c=200  !  number  of 

T  au^ 

VALUES 

20 

Dtc=.01  1  INCREMENT 

IN  Tau^ 

30 

N=1024  !  SIZE  OF  FFT;  N 

>  Nt  c 

REQUIRED 

40 

DOUBLE  Htc ,N,N4,N2,H£  >  INTEGERS 

50 

N4  =  N/4 

60 

N2=N^2 

70 

REDIM  Co£<0; N4) , X<0; H-1) , Y(0:N-1 > 

80 

DIN  Co£<256) , X< 1023) , Y< 1023), fi< 200) 

90 

T=2. *PI/N 

100 

FOR  N£=0  TO  N4 

110 

Cos(:N£)=C0S<T*Ns)  !  QUfiRTER-COSINE 

TABLE 

IN  Cost*) 

120 

NEXT  Ns 

130 

GINIT  200/260 

140 

PLOTTER  IS  505, "HPGL" 

150 

PRINTER  IS  505 

160 

LIMIT  PLOTTER  505,0,200,0,260 

170 

VIENPORT  22,85, 19, 122 

180 

WINDOW  0,N2,-5, 1 

190 

PRINT  "VS5" 

200 

GRID  H/10, 1 

210 

PRINT  "VS36" 

220 

ASSIGN  #1  TO  "fllTl" 

230 

READ  IH;A<*) 

240 

NAT  X=<0. ) 

250 

MAT  Y=<:0.) 

260 

X(:0)  =  .  5*A<0) 

270 

FOR  H£=1  to  Ntc-1 

280 

X<N£)=A<N£) 

2  :i0 

NEXT  Ns 

300 

XCNlc  >  =  . 5*A<Nlc  ) 
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310 

320 

330 

340 

350 

360 

370 

380 

390 

400 

410 

420 

430 

440 

450 

460 

470 

480 

490 

500 

510 

520 

530 

540 

550 

560 

570 
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MAT  X=X*<Dtc*4.) 

CALL  Fft  14kH,Co£C*>,X<:*>,Y<*)) 
FOR  Hs=0  TO  H2 
Ar  =  X<:H£> 

IF  Ar>0.  THEN  380 

PENUP 

GOTO  390 

PLOT  Hs.LGTCAr) 

NEXT  As 
PEHUP 

ASSIGN  #1  TO  "AOT0" 

READ  #1;A(*) 

MAT  X=<0. ) 

MAT  Y=<0. ) 

X<0>=. 5*A<0) 

FOR  Ns=l  TO  Htc-1 
X(:N£)=A<Ns) 

NEXT  Ns 

X<Ntc)=.5*A<Ntc> 

MAT  X  =  X*<Dtc*4.  ) 

CALL  Fft  14<N,Co£<*),X<*),Y<:*)) 

FOR  N£=0  to  N2 

Ar=X<Ns) 

IF  Ar>0.  THEN  570 

PENUP 

goto  580 

PLOT  Ns,LGT<Ar) 

NEXT  Ns 
PENUP 

ASSIGN  #1  TO  "AOTl" 

READ  #1;A<*) 

NAT  X=<0. ) 

MAT  Y=<0. ) 

X<0)=. 5*A<0) 

FOR  Ns=l  TO  Ntc-1 
X<N£)=A(Ns) 

NEXT  Ns 

X<Ntc )=. 5*A<Ntc ) 

MAT  X  =  X*<Dtc*4.  ) 

CALL  Fftl4<N,Cos<*),X<*),Y<*)) 

FOR  Ns^O  TO  N2 

Ar=X<Ns) 

IF  Ar>0.  THEN  760 

PENUP 

GOTO  770 

PLOT  Ns,LGT(Ar) 

NEXT  Ns 
PENUP 

ASSIGN  #1  TO  "A0T2" 

READ  #1;A<*) 

MAT  A=A/(A<0)) 

MAT  X=<0.  ) 

NAT  Y=<0.> 
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840  X<0>=.5*fl<0> 

850  FOR  Ns=l  TO  Htc-1 

860  XcM£>=FI':Hs> 

870  NEXT  Us 

880  XCHtc )=. 5*flCNlc > 

890  MAT  X=X*<Dlc*4.) 

900  CALL  Fftl4<H,Cos<*),X(*),Y<*)> 

910  FOR  Ns=0  TO  N2 

920  Ar=X<N5) 

930  IF  Ar>0.  THEN  960 

940  PENUP 

950  GOTO  970 

960  PLOT  Ns,LGT(Rr) 

97<s  NEXT  Ns 

980  PENUP 

990  PAUSE 

1000  END 

1010  ! 

1020  SUE  Fft.  14<D0UBLE  N.REAL  Cos  <  *  )  ,  X  (  *  )  ,  Y  (  *  )  )  !  N<  =2'- 1  4=  1 6384  ;  0  SUES 

1030  DOUBLE  Log2n,Nl,N2,N3,N4,  J,K  !  INTEGERS  <  2'-31  =  2,147,483,648 

1040  DOUBLE  II,  1 2, 13, 14, 15, 16, 17, 18, 19, 110, Ill, 112, 113, 114, L<0:13> 

1050  IF  N=1  THEN  SUEEXIT 

1060  IF  N>2  THEN  1140 

1070  A=X<0)+X<1) 

1083  X< 1 )=X<0)-X< 1 ) 

1090  X<0)=A 

1100  A=Y<0)+Y<1> 

1110  Y< 1 )=Y<0>-Y<  1  ) 

1120  Y<0)=A 

1130  SUBEXIT 

1140  A=L0G<N)xL0G<2. ) 

1150  Log2n=R 

1160  IF  ABS<A-Log2n>< 1 . E-8  THEN  1190 

1170  PRINT  "N  ="5N5"IS  HOT  A  POWER  OF  2;  DISALLOWED." 

1180  PAUSE 

1190  Nl=N/4 

1200  N2=N1+1 

1210  N3=N2+1 

1220  N4=N3+N1 

1230  FOR  11=1  TO  Log2n 

1240  12  =  2''<Log2n-1 1  ) 

1250  13=2*12 

1260  I4=N/I3 

1270  FOR  15=1  TO  12 

1280  I6=< 15-1 )*I4+1 

1290  IF  I6<=N2  THEN  1330 

1300  A1=-Co£<N4-I6-1 ) 

1310  A2=-Co£< I6-N1-1 ) 

1320  GOTO  1350 

1330  Al=Cos<I6-l) 

1340  R2  =  -Co£<N3-I6-1  ) 

1350  FOR  17=0  TO  N-I3  STEP  13 

1360  18=17+15-1 

1370  19=18+12 

1380  T1=X«:I6) 

1390  T2=X<I9) 
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T3  =  Y<:  I8> 

T4  =  Y<:  19) 

Fl3  =  Tl-T2 

fi4=T3-T4 

X( IS)=T1+T2 

Y( I8)=T3+T4 

X( I9)=fll*fl3-fl2*fi4 

YC  I9)=fll*fl4  +  n2*fi3 

NEXT  17 

NEXT  T5 

NEXT  I  I 

I  1  =Log2ri+  1 

FOR  12=1  TO  14 

L< 12-1  )  =  1 

IF  !2>Log2ri  THEN  1560 
L<  12-1  )  =  2''<  11-12) 

NEXT  12 
K  =  0 

FOR  11=1  TO  L<13) 

FOR  12=11  TO  La2)  STEP  L<13) 
FOR  13=12  TO  LCll)  STEP  L<12) 
FOR  14=13  TO  L<10)  STEP  L<11) 
FOR  15=14  TO  L<9)  STEP  L<10) 
FOR  16=15  TC  L<3)  STEP  L<:9) 
FOR  17=16  TO  L<7)  STEP  L<8) 
FOR  18=17  TO  L<6)  STEP  L<7) 
FOR  19=18  TO  L<5)  STEP  L<6) 
FOR  110=19  TO  L<4)  STEP  LC5) 


FOR 

111  =  110 

TO 

L<3) 

STEP 

L<4> 

FOR 

112=111 

TO 

L<2) 

STEP 

L<3) 

FOR 

113=112 

TO 

L<  1  ) 

STEP 

L<2) 

FOR 

114=113 

TO 

L<0) 

STEP 

L<1) 

J=1 14-1 

IF 

K>J  THEN 

1800 

fl=X(K) 
X<K)=X< J) 
XC J)=fl 
fl=Y<K) 
Y<K)=Y< J) 
Y< J)=fi 
K.  =  K+1 
NEXT  114 
NEXT  113 
NEXT  112 
NEXT  Ill 
NEXT  110 
NEXT  19 
NEXT  18 
NEXT  17 
NEXT  16 
NEXT  15 
NEXT  14 
NEXT  13 
NEXT  12 
NEXT  I  1 
SUBEND 
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APPENDIX  A. 4  -  EVALUATION  OF  WAVENUMBER  INTENSITY 
SPECTRUM  FOR  ZERO  DELAY  (f,f'  =  0) 


A  program  for  the  numerical  evaluation  of  the  zeroth-order 
Bessel  transform  of  covariance  My(AR,0)  -  is  contained  in 

this  appendix.  Inputs  required  of  the  user  are  listed  in  lines 
10  -  40  and  are  coupled  to  appendix  A. 2,  where  the  data  input, 
AORO  or  AORl  or  AOR2 ,  was  generated.  The  numerical  Bessel 
transform  is  accomplished  by  means  of  Simpson's  rule  with  end 
correction  [11;  pages  414  -  418],  and  is  exceedingly  accurate  for 
the  small  increment,  .005,  in  AR  employed  in  line  30. 
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Dkc  =  .  4 
Hkc=200 
Dr-c  = .  005 
Mrc=900 

DOUBLE  Mrc,Mkc,I,Ns 
REDIM  C<0:Hrc) 

REDIM  Hi <0:Hkc>,W0<0:Nk 
DIM  C<900) , Wi <200) , W0<2 
ASSIGN  #1  TO  "filRl" 

READ  #1;C<*) 

FOR  1=0  TO  Nkc 
Kc=I*Dkc 
T  =  Kc  *Drc 


!  INCREMENT  IN  k^ 

!  HUMBER  OF  k-^  VALUES 
!  INCREMENT  IN  DelR^ 

!  NUMBER  OF  DelR^  VALUES 
!  INTEGERS 

>,Wl<0:Nkc),W2<0;Nkc> 

0),W1<200),W2<200) 


!  k'' 


Se=So=0. 

FOR  Ns=l  TO  Hrc-1  STEP  2 
So  =  So  +  Ns*FNJo<T*Ns)*C  <N£) 
NEXT  Ns 

FOR  Ns=2  TO  Nrc-2  STEP  2 
Se=Se+Ns*FHJo(T*Ns)*C<Ns) 
NEXT  Ns 


Wi < I )=C<0)+16. *So+14. *Se 
NEXT  I 

MAT  Wi  =Wi  *<Drc*Lrc*2.  *PI/'15.  ) 
ASSIGN  #1  TO  "AOR0" 

READ  tHjC<*) 

FOR  1=0  TO  Nkc 
Kc  =  I *Dk  c 
T=Kc*Drc 


Se=So=0 . 


FOR  Ns=l  TO  Nrc-1  STEP  2 


So  =  So  +  NS'»FNJo<T*Ns)*C<Ns) 


NEXT  Ns 

FOR  Ns=2  TO  Nrc-2  STEP  2 
Se=Se+Ns*FNJo<T*Ns)*C<Ns) 


NEXT  Ns 

W0< I )=C<0>+16. *So+14. ♦Se 


NEXT  I 

NAT  W0  =  W0*<Drc*Drc*2.  *PI-'15.  ) 
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390  fiSSIGU  #1  TO  "RORl  ' 

400  READ 

410  FOR  1=0  TO  Hkc 

420  Kc=I*Dkc 

430  T=Kc*Drc 

440  Se=So=0. 

450  FOR  Ms=l  TO  Hrc-1  STEP  2 

460  So  =  So  +  N£*FMJo<  T*H£ )*C<Hs) 

470  NEXT  Ns 

480  FOR  Ns=2  TO  Nrc-2  STEP  2 

490  Se=Se+N£*FN Jo<T*Ns)*C(Ns) 

500  NEXT  Ns 

510  HI < I >=C<0>+16. *Go+14. »Se 

520  NEXT  I 

530  MAT  Ul=Nl»CD('C*Drc*2.*PI/15.  ) 

540  ASSIGN  #1  TO  "A0R2" 

550  READ  #1;C<:*) 

560  ASSIGN  #1  TO  * 

570  FOR  1=0  TO  Nkc 

580  Kc=I*Dkc 

590  T=Kc*Drc 

600  Se=So=0. 

610  FOR  Ns=l  TO  Nrc-1  STEP  2 

620  So=So+Ns*FNJo<T*Ns)*C<Ns) 

630  NEXT  Ns 

640  FOR  Ns=2  TO  Nrc-2  STEP  2 

650  S€-  =  Se  +  N£^r-  ■!Jo<:T*Ns)*C<Ns> 

660  NEXT  Ns 

670  W2<n=C<e:)  +  16.*So+14.«Se 

680  NEXT  I 

690  MAT  W2=W2*<Drc*Drc*2. *PI/15. > 

700  GINIT  200/260 

710  PLOTTER  IS  505,"HPGL" 

720  PRINTER  IS  505 

730  LIMIT  PLOTTER  505,0,200,0,260 

740  VIEWPORT  22,85,19,122 

750  WINDOW  0, Nkc , -9, 1 

760  PRINT  "VS5" 

770  GRID  25,1 

780  PRINT  -'0336" 

790  FOR  1=0  TO  Nkc 

800  W  =  Wi<n 

810  IF  W>0.  THEN  840 

820  PENUP 

830  GOTO  850 

840  PLOT  I,LGT<W) 

850  NEXT  I 

860  PENUP 
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870  FOR  1=0  TO  Hkc 

880  W=W0CI> 

S90  IF  IJ>0.  THEN  920 

900  PENUP 

910  GOTO  930 

920  PLOT  I,LGT<W) 

930  NEXT  I 

940  PENUP 

950  FOR  1=0  TO  Nkc 

960  W=W1(I) 

970  IF  W>0.  THEN  1000 

980  PENUP 

990  GOTO  1010 

>.000  PLOT  I,LGT<W> 

1010  NEXT  I 

1020  PENUP 

1030  FOR  1=0  TO  Nkc 

1040  W  =  W2<:i) 

1050  IF  W>0.  THEN  1080 

1060  PENUP 

1070  GOTO  1090 

1080  PLOT  I,LGT<W) 

1090  NEXT  I 

1100  PENUP 

1110  PAUSE 

1120  PRINTER  IS  CRT 

1130  PLOTTER  505  IS  TERMINATED 

1140  END 

1150  ! 

1160  DEF  FNJo<X)  !  Jo<X)  FOR  ALL  X 

1170  Y=RBS<X) 

1180  IF  Y>8.  THEN  1280 

1190  T=Y*Y  !  HART,  #5845 

1200  P  =  227 1  490439. 5536033-T*< 55 13584. 5647707522-T*5292. 6171303845574:) 

1210  P  =  2334  4  891  71  877869. 7-T*  <  47765559442673 . 588-T*  (  462 1  7222503 1 . 71803-T*p:) 

1220  P= 1859623 17621 897804. -T*<44 14582939 181 5982. -T*P> 

1230  0  =  204251483. 52 134357  +  T*< 494030.  7949 1 8 1  3972  +  T* ( 884 . 72036756  1  75504  +  T  > !) 

1240  Q  =  254475001 3658996. 8+T»<  15015462449769. 752+T*< 64398674535.  133256  +  T*G!  :> 

1250  0=185962317621897733. +T*Q 

1260  Jo=P/Q  ' 

1270  RETURN  Jo 

1280  2  =  8. /Y  !  HART,  #6546  8.  6946 

1290  T=2»Z 

1300  Pri=2204. 5010439651804+T*< 128. 67758574871419+T*. 90047934748028803) 

1310  Pn=8554. 82254 150666 17+ T*< 8894. 4375329606 194 +T*Pn> 

1320  Pd=22 14. 048851 91471 04+T»< 130. 88490049992388+T) 

1330  Pd=8554. 82254 15066628+ T*< 8903. 8361 4 17095954 +T*Pd> 

1340  Qn=13. 99097686596O680+T»< 1 .0497327982345548+T*. 0G935259532940319> 

1  350  Qn  =  -37. 510534954  9571 12-T+C 46. 0938268 14625 175+T*0ri) 

1360  Od=921 . 5669755265309O+T*<74. 42838974 141 1 179+T> 

1370  Od=2400. 6742371 1 72675+T* <2971 . 9837452084920+T+Qd > 

1380  T=Y-. 73539816339744828 

1390  Jo=.  28209479 177387S20*SOR<2)»< COS <T>*Pn-'Pd-SIN<T:)*2*eri.-0d;) 

1400  RETURN  Jo 

1410  FNEND 
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APPENDIX  A. 5  -  EVALUATION  OF  PHASE  MODULATION  INTENSITY  SPECTRUM 

The  normalized  covariance  function  for  phase  modulation  is 
given  by  (2.50)  in  the  main  text,  namely 

k  (C.)  =  expf-  r;  [1  -  expi-:,)]  -  A[2  -  p(2.)]  + 

°  LAP  (A. 5-1) 

^2 

+  2A  [1  -  p(C)]  exp(-p^/A)  +  A  p(J.)  exp(-  [1  -  exp(-bC)])] 

for  C  i  0,  where  C  is  the  time  delay  and  p(C)  is  the  temporal 

2  2 

normalized  covariance  of  the  field  process.  Also  Pp  =  PpQ- 
Since  (A. 5-1)  involves  an  exponential  of  an  exponential  of  an 
exponential,  and  because  a  wide  range  of  parameter  values  are  of 
interest,  care  must  be  taken  in  numerical  evaluation  of  this 
covariance  and  its  transform. 

Observe  first  that 

k^(0)  =  1  since  p{0)  =  1  .  (A, 5-2) 

Also,  as  delay  C  +®,  then  p  -»  0,  giving 

k^C®)  =  exp[-r^  /jp  -  2A  +  2A  exp(-pp/A]]  ^  0  .  (A. 5-3) 

The  spectrum  of  interest  is  given  by 


Wq(w)  =  4  J  dC  cos(wC,)  k^(C)  for  w  i  0  ;  w  =  2nf  .  (A. 5-4) 

0 

The  nonzero  value  of  (A. 5-3)  at  C  -  *  leads  to  an  impulse  in 
spectrum  W^(a))  at  u  =  0.  This  limiting  value,  k^(®),  must  be 
subtracted  from  covariance  (A. 5-1)  prior  to  the  numerical 
Fourier  transform  indicated  by  (A.5-4). 
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2 

For  /Jp  <<1/  the  term 

expj-  pp  [1  -  exp{-C.)])  (A.5-5) 

approaches  its  limiting  value  at  C,  =  +“  as  follows: 

exp(-  pp  [1  -  exp(-C)])  -  exp(-  Pp)  = 

=  exp(-  pp]  [exp(r^  Pp  exp(-C))  -  1  = 

~  exp(-  Pp]  Pp  exp(-:.)  .  (A. 5-6) 

This  is  a  fairly  rapid  decay  with  Z,  and  will  not  lead  to 

2 

numerical  difficulty  when  Pp  <<  1. 

For  large  bpp/A,  the  term 

exp[ - 1  ■  exp(-bC.)]]  (A. 5-7) 

is  very  sharp  near  C  =  0;  in  fact,  it  is  given  approximately  by 

^2 

exp|-  bcj  for  Z,  near  0  .  (A. 5-8) 

Therefore,  we  define  the  sharp  component  of  covariance  as 

^2 

kg(C)  =  exp|-  A  +  A  exp|-  bf,j  I  -  exp(-A)  for  all  Z,  .  (A. 5-9) 
Then 

k^(0)  =  1  -  exp(-A)  ,  k  (®)  =  0  .  (A. 5-10) 

s  & 
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Now  we  let 


=  +  kg(2:)  ,  (A. 5-11) 

where  k^(C,)  is  a  flat  function  near  2,  *  0.  Then  we  can  express 
the  desired  difference  as 

ko(C)  -  k^(«o)  «  [k^(J.)  -  k^C*)]  +  kg(C)  « 

2  k,^(C)  +  kg(2.)  /  (A. 5-12) 

where  functions  k,  (J.)  and  k  {^)  both  decay  to  0  at  ^  We  now 

A  9 

employ  two  separate  FFTs  on  each  of  the  functions  in  (A. 5-12). 

The  sharp  component,  k  (2.),  must  be  sampled  with  a  very  small 

2 

increment,  A2,,  when  bpp/A  is  large.  On  the  other  hand,  the  flat 
component 

k^(C)  '  kj(C)  -  k^C)  {A.5-13) 

2 

can  be  sampled  in  a  coarser  fashion.  Finally,  if  bpp/A  is 
moderate,  we  work  directly  with  kQ(C)  -  k^(«»)  without  breaking 
it  into  any  components. 

Two  programs  are  furnished  in  this  appendix,  one  for  moderate 

2  2 
bpp/A,  and  the  other  for  the  flat  component  (A.5-13)  when  bpp/A 

is  large.  For  sake  of  brevity,  the  Fourier  transform  of  the 

sharp  component  (A. 5-9)  is  straightforward  and  is  not  presented. 

The  particular  covariance  p(C)  adopted  is  triangular, 

p(C)  *  1  -  for  Id  <  C-,  »  0  otherwise  ,  (A. 5-14) 

but  can  easily  be  replaced.  The  parameter  is  the  cutoff  value 
of  covariance  p{d* 
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The  number  of  samples,  N,  taken  of  the  covariance,  in  order 
to  perform  the  FFT  of  (A.5-4),  is  rather  large,  so  as  to 
guarantee  a  very  small  value  of  truncation  error  at  the  upper  end 
of  the  integral,  despite  the  small  increment  AC*  In  order  to 
keep  the  FFT  size,  Mf,  at  reasonable  values,  the  data  sequence  is 
collapsed  without  any  loss  of  accuracy  according  to  the  method 
given  in  [8;  pages  7-8]  and  [9;  pages  13  -  16].  The 
trapezoidal  rule  is  used  to  approximate  the  integral  in  (A.5-4), 
for  reasons  given  in  [8;  appendix  A]. 
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SPELTRUtl  FOR  PHASE  MODULATION  -  MODERATE 


M  u  R  =  1  .  ! 

Gr=.001  I 

Bs=l.  ! 

A=.  2  ! 

2c=2.*PI  I 

De1z=.005  ! 

H=60000  ! 

Mf= 16384  ! 

DOUBLE  N,Mf,Ms,Ms  ! 

DIM  X';  1630^’  > ,  Y<  16384)  ,  Cos 
RED  I M  X  C  0 ; Mf - 1 ) , Y  <  0 ; Mf- 1 ) 
MAT  X=(0.)  . 

MAT  Y=<0.  ) 

T  =  2.'»PI/Mf 


MUsubP 

Gamma" 

b 

A 

Rho<Z)  =  0  for  lZ|>Zc; 
Zeta  increment 
Maximum  number  of  samp 
Size  of  FFT 
INTEGERS 
4096) 

C:os<0;Mfx4) 


Z  =  z  e  t  a 

e  s  o  f  k  o  ( z  e  t  a ) 


FOR  M£=0  to  Mf/4 

Cos<:M£)=C0S(T*Ms)  !  QUARTER-COSINE  TABLE 

NEXT  Ms 

T  a=Gp»MuR*Mup 

IF  A=0.  THEN  220 

Tb  =  Mup*Mup'"A 

Tc=2. *R*FNExp(Tb) 

Ki nf=FNExp(Ta+2. *A-Tc )  !  CORRELATION  AT  INFINITY 

COM  A,Bs,Zc,Ta,Tb,Tc,Kinf 
T=1 . -Ki nf 
PRINT  0,T 

X<0)=T*.5  !  TRAPEZOIDAL  RULE 


FOR  N£=1  to  N 

Corr=FNKo < N£*De 1 z )  !  CORRELATION  koCzeta) 

IF  Ns<6  THEN  PRINT  Ns.Corr 
IF  ABS<Corr)< 1 . E-30  THEN  350 
Ms=Ns  MODULO  Mf  !  COLLAPSING 

X  <  Ms ) =X  <  Ms ) +Cor r 


NEXT  Ns 

PRINT  "Final  value  of  Corr  =";Corr;"  Ns  =";Ns 
MAT  X=X»(Del z*4. ) 

CALL  Fft  14<Mf  ,Cos<*),X<*),  Y<^«-)) 
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380  GINIT 

390  PLOTTER  IS  "GRAPHICS" 

400  GRAPHICS  OrJ 

410  NIMHOH  -2, 2, -60,0 

420  LINE  TYPE  3 

430  GRID  1,10 

440  LINE  TYPE  1 

450  Delf=l./<Mf*Delz) 

460  FOR  Mi=l  TO  Mf-"2 

470  F=Ms»De1f  !  FREQUENCY 

480  PLOT  LGT<F) , 10. *LGTCX<Hs) ) 

490  NEXT  Ms 

500  PENUP 

510  PAUSE 

520  END 

530  ! 

540  DEF  FNExp<Xminus>  !  EXP<-X>  WITHOUT  UNDERFLOW 

550  IF  Xfiii  r.us>?08. 3  THEN  RETURN  0. 

560  RETURN  EXP<-Xminu£) 

570  FNEND 

580  ! 

590  DEF  FNKoCZeta)  !  CORRELATION  koCzeta) 

600  CON  A, Bs , 2c , Ta, Tb, Tc , Ki nf 

610  Rho*NAX(0. , 1 . -2et a/2c >  !  TRIANGULAR  RHO 

620  El=Ta*( 1 . -FNExpCZeta) ) 

630  E2=Tb*< 1 . -FHExp<Bs*2eta>) 

640  E3=A*Rhc*FNExp<E2) 

650  RETURN  FHExp < E 1 +A* < 2 . -Rho ) -Tc * < 1 . -Rho > -E3 > -K i nf 
660  FNEND 

670  ! 

680  SUE  Fft  14<D0UBLE  H,RERL  Cos(*> ,  X<») ,  Y<»> )  !  N<  =2'- 1 4=  1 6384  ;  O  SUE 
690  !  SEE  APPENDIX  A. 3 
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SPECTRUM  FOR  PHASE  MODULATION  -  FLAT  COMPONENT 


N  u  p  =  1 . 

Gp= . 00 1 
Bs=l  . 

A  =  0. 

2c=2. *PI 
Delz=.005 
N=60000 
Mf= 16384 

DOUBLE  N,Mf,Ms,Ns 
DIN  X<16384),Y<16384),Cos<4096) 

REDIN  X<0: Mf-1 >, Y<0! Mf-l ),Cos(0;Mrx4) 

NAT  X=<0. ) 

MAT  Y=(0. ) 

•'  to  •  •  A  '  I  •  t 

FOR  M£=0  TO  Mfx4 

Cos<Ms)=COS<T*H£>  !  QUARTER-COSINE  TABLE 

NEXT  Ms 


MUsubP 

Gamma'' 

b 

A 

Rho<2)  =  0  for  |Z|>2c;  2=zeta 
Zeta  increment 

Maximufri  number  of  samples  of  kl(zeta) 

Size  of  FFT 

INTEGERS 


T  a=Gp*Mup»Mup 
IF  A=0.  THEN  220 
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210  Tb  =  tUip*Mup 'FI 

220  Tc=2. *ft*FHExp(Tb) 

230  Tb=5.E55 

240  Ki  nf  =  FHExp<T=i+2.  *fi-Tc  >  !  CORFELFlT  I  OH  AT  INFINITY 

250  Ea=FNExp<R> 

2e0  Tbb=Tb*Bs 

270  COM  A,Bs,2c,Ta,Tb,Tc,Kinf,Ea,Tbb 

280  T=1 . -Ki nf-( 1 . -Ea>  !  SUBTRACT  SHARP  COMPONENT 

290  PRINT  0,T 

300  X'.'0>=T*.5  !  TRAPEZOIDAL  RULE 

310  FOR  Ns=l  TO  N 

320  Corr=FNKl CNs^Delz)  !  CORRELATION  LlCzet a> 

330  IF  Ns<6  THEN  PRINT  Ns, Corn 

340  IF  ABSCCorrX  1  .  E-30  THEN  380 

350  M£=Ns  MODULO  Mf  !  COLLAPSING 

360  X(Ni)=X':Ms>+Corr' 

370  NEXT  Ns 

380  PRINT  "Final  ualue  of  Corn  =";Corr;"  Ns  =";Ns 
390  MAT  X=X*<De 1 2*4. > 

400  CALL  Ffl  14CMf,Cos<*),X<*),Y(*'-«) 

410  G I H I  T 

420  PLOTTER  IS  "GRAPHICS" 

430  GRAPHICS  ON 

440  WINDOW  -2, 2, -60,0 

450  LINE  TYPE  3 

460  GRID  1,10 

470  LINE  TYPE  1 

480  De 1 f=l . /<Nf*Del z) 

490  FOR  Ms=l  Tn  Mf/2 

500  F=Ms*Delf  !  FREQUENCY 

510  T=X<Ms) 

520  IF  T>0.  THEN  550 

530  PENUP 

540  GOTO  560 

550  PLOT  LGT<F) , 10. *LGT<T) 

560  NEXT  Ms 

570  PENUP 

580  PAUSE 

590  END 

600  ! 

610  DEF  FNExpOitrii  nus)  !  EXP<-X>  WITHOUT  UNDERFLOW 

620  IF  Xfni  nus>708. 3  THEN  RETURN  0. 

630  RETURN  EXPC-Xminus) 

640  FNEND 

650  ! 

660  DEF  FNKlCZeta)  !  CORRELATION  klCzeta) 

670  COM  A,Bs,2c,Ta,Tb,Tc,Kinf,Ea,Tbb 

680  Rho  =  MAX<:0.  ,  1  . -2eta/2c )  !  TRIANGULAR  RHO 

690  El=Ta*( 1 . -FNExpCZela) ) 

700  E2=Tb*< 1 . -FNExp<Bs*2eta) ) 

710  E3=A»Rho*FNExp<E2) 

720  E4=FNExp<Tbb*Zeta> 

730  Shar p=FNExp C A* < 1 . -E4 ) ) -Ea  !  ks<Z£la) 

740  RETURN  FNExp < E 1 +A* < 2 . -Rho) -Tc * ( 1 . -Rho > -E3 > -K i nf -Shar p 
750  FNEND 

760  ! 

770  SUE  Fftl4<DOUBLE  N.REAL  Cos < * ) , X < * ) , Y ( * ) >  !  N< =2^ 1 4= 1 6384 ;  0  SUBS 

780  !  SEE  APPENDIX  A. 3 
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APPENDIX  A. 6  -  EVALUATION  OF  FREQUENCY  MODULATION 
INTENSITY  SPECTRUM 

The  normalized  covariance  function  for  frequency  modulation 
is  given  by  (2.48)  in  the  main  text,  namely 

k^(C)  =  exp[-  [exp(-2.)  +  C  -  1]  -  A[2  -  p(C)]  + 

^2 

+  A  p{C)  exp[-  [exp(-bC)  +  bC  -  1])]  for  C  ^  0  ,  (A.6-i) 

where  L  is  the  time  delay  and  p(C)  is  the  temporal  normalized 

2  2 

covariance  of  the  field  process.  Also,  Pp  =  pp^  and 
b  =  Since  (A. 6-1)  involves  an  exponential  of  an 

exponential  of  an  exponential,  and  because  a  wide  range  of 
parameter  values  are  of  interest,  care  must  be  taken  in  numerical 
evaluation  of  this  covariance  and  its  transform. 

Observe  that 

^o^^^  ~  ^  '  because  p(0)  =  1  .  (A. 6-2) 

Also,  as  delay  Z,  -*  +’*>,  then  p  -»  0,  giving 

k^CC)  ~  exp(-  p^  (C  -  1)  -  2a]  h  kj(C)  for  C  >  0  .  (A. 6-3) 

2 

This  term,  k^{C)/  decays  slowly  with  C  if  Pp  <<  1  . 

The  spectrum  of  interest  is  given  by 

00 

Wq(w)  *  4  J  dC,  cos(a)C)  w  2  0  ;  w  =  2nf  .  (A. 6-4) 

0 

The  spectrum  corresponding  to  the  limiting  component,  kj^{C,)  in 
(A. 6-3),  is  directly  available  in  closed  form  as 
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Wj(a>)  =  4  J  dC  cos(u)C)  kj(C.)  *= 
0 


*  exp 


4  r'  p 


2 

F _ 

+  fa) 


2  * 


(A. 6-5) 


2 

If  /Dp  <<  1,  this  latter  quantity  is  large  and  very  sharply 

peaked  at  (o  =  0;  hence,  this  term  has  been  subtracted  out  and 

2 

handled  separately  when  Pp  <<  1*  The  residual  covariance, 
k^(C,)  -  kj(C.),  then  decays  very  rapidly  with  C.  and  is  easily 
handled  directly  by  means  of  an  FFT.  This  breakdown  is  not 
necessary  when  Pp  ~  1  and  is  avoided,  then,  by  handling 
kQ(C,)  directly  in  one  FFT. 

For  /Jp/A  >>  1,  the  term 

f  Pp  ^ 

exp - =■  [exp(-bC)  +  bC,  -  1]  (A. 6-6) 

^  Ab*^  ^ 

inside  the  exponential  in  (A. 6-1)  behaves  like 

^2 

exp( - I  ^  near  C  *  0  ,  (A. 6-7) 

where  its  major  sharp  contribution  arises.  For  example,  if 
M  *  50,  A  =  1,  then  increment  AC  ■  .005  leads  to  values  for 

T 

(A. 6-7)  of 

exp( -0.156  n^)  at  C  *  n  AC  ,  (A. 6-8) 

which  is  adequately  sampled  in  order  to  track  its  dominant 
contribution;  the  actual  sequence  of  values  is  1,  .856,  .536, 
.247,  .083.  For  smaller  Pp/A,  this  sampling  interval  is  more 
than  adequate. 
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Two  programs  are  furnished  in  this  appendix,  one  each  for  the 

2 

cases  of  large  and  small  Pp.  The  particular  covariance  p(C) 
adopted  is  triangular, 

p(C)  =  1  -  for  Id  <  /  0  otherwise  ,  (A. 6-9) 

^c 

but  can  easily  be  replaced.  The  parameter  is  the  cutoff 
covariance  value  and  is  specified  numerically  in  the  examples  in 
figures  3.9a  through  3.10b. 


1  0 

!  SPECTRUM  FOR  FREQUENCY  MODULfiTION  -  LfiRGE  Gp*Muf--2 

20 

Muf=50.  ! 

MUsubF 

30 

G(j=.0Ol  ! 

G  atn  rn  a  ^ 

40 

Bs=l.  ! 

b 

50 

Ft=.  2  ! 

A 

60 

2c=2.*PI  ! 

Rho<2)  =  0  for  |2|>2c;  Z  =  zet.a 

70 

Il€l2=.005  ! 

2eta  increme-nt 

80 

H=6O000  ! 

N ax i m u m  number  of  s am p 1 e s  o f  k  o < z e t  a ) 

90 

Mf= 16384  ! 

Size  of  FFT 

100 

nOUBLE  N,Mf,M£,Ns  ! 

INTEGERS 

110 

DIM  X<  163ci4) ,  Y(  16384) ,  Cos <4096) 

120 

REDIM  X<0: Mf -1 ) , Y<0: Mf-1 ) , 

Cos<0!Mf/'4) 

130 

MAT  X=<0. ) 

140 

MAT  Y=<0. ) 

150 

T=2. *PI/Mf 

160 

FOR  Ms  =  0  TO  rif/'4 

170 

Cos<Ms)=C0S<T^Ms)  ! 

QUARTER-COSINE  TABLE 

180 

NEXT  Ms 

190 

T  a=Gp-5‘Muf  ♦Muf 

200 

IF  A=0.  THEN  220 

210 

Tb  =  Muf*Muf/'<A*Bs*Bs) 

220 

Tc=FNExp<2.  *A-T%)*Tck 

230 

T  d  =  T  a*T  a 

240 

COM  H , Bs , 2c , Ta, Tb 

250 

X<0)=.5  1 

TRAPEZOIDAL  RULE 

260 

FOR  Ns=l  TO  N 

270 

Corr=FHKo<Ns*Delz)  ! 

CORRELATION  ko<zeta) 

280 

IF  Corr<l.E-20  THEN  320 

290 

Ms=Ns  MODULO  Mf  ! 

COLLAPSING 

300 

X<Ms)=X<Ms)+Corr 

310 

NEXT  Ns 

320 

PRINT  "Final  ualue  of  Corn 

=  '’;Corr;"  Ns  =";Ns 

330 

NAT  X=X*<Delz) 

340 

CALL  Fft 14<Nf ,Cos<»),X<*), 

Y  <  *  )  ) 
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350  GINIT 

3i;0  PLOTTER  IS  "GRlTPHICS" 

370  GRFiPHICS  OH 

380  WINDOW  -4,2,-70,30 

390  LINE  TYPE  3 

400  GRID  1,10 

410  LINE  TYPE  1 

420  De  1  ^=1  . /01P*Dt  1  z  ) 

430  FOR  Ms=l  TO  Mf.'2 

440  F=ns*Delf  !  FREQUENCY 

450  T  =  X<ri5) 

460  IF  T>0.  THEN  490 

470  PENUP 

480  GOTO  500 

490  PLOT  LGTCF), 10.*LGTCT) 

500  NEXT  Ms 

510  PENUP 

520  Fldd  =  X<:0)-Tc/Td  !  ORIGIN  CORRECTION 

530  F=l.E-4 

540  FOR  Ns=l  TO  100 

550  W  =  2.itPI*F 

560  Wl=Tc/'<Td  +  W*W) 

570  T=Wl+fldd 

580  IF  T>0.  THEN  610 

590  PENUP 

600  GOTO  620 

610  PLOT  LGT(F>, 10.*LGT<Wl+fldd) 

620  F=F*1.1  !  FREQUENCY 

630  IF  F>Delf  THEN  650 

640  NEXT  Ns 

650  PENUP 

660  PAUSE 

670  END 

680  ! 

690  DEF  FNExp<Xmi nus)  !  EXP(-X)  WITHOUT  UNDERFLOW 

700  IF  Xfni  nus>708. 3  THEN  RETURN  0. 

710  RETURN  EXPC-Xminus) 

720  FNEND 

730  ! 

740  DEF  FNKoCZela)  !  CORRELATION  ko(zeta) 

750  COM  A, Bs, 2c , Ta, Tb 

760  El=FNExp<2eta)+2eta-l . 

770  T=Bs*Zeta 

780  E2=FNExp<T)+T-l . 

790  Rho  =  MAX(:0.  ,  l-2etax2c  )  !  TRIANGULAR  RHO 

800  T  =  Ta»E  1  +A*<:2.  -  Rho)-A»Rho*FNExp<Tb*E2) 

810  RETURN  FNExpCT) 

820  FNEND 

830  ! 

840  SUB  Fft 14 (DOUBLE  N,REAL  Cos < » > , X( * > , Y < * > >  !  N< =2  1 4= 1 6384 ;  O  SUBS 

850  !  SEE  APPENDIX  A, 3 
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SPECTRUM  FOR  FREQUENCY  MODULFiT  I  ON  -  SMRLL  Gp*Nuf''2 
Muf=l.  !  MUiubF 

Gfj=.001  i  Gatnma^ 

Bs=l.  !  b 

B=.2  !  fl 

2c=2.*PI  !  Rho<2)  =  0  for  | 2 | >2c ; 

Del  2=.  005  !  Zeta  incretnent 

N=  10000  !  Maxi  mum  riumber  of  sampl 

Hf=S192  !  Size  of  FFT 

INTEGERS 


2  I  >2c 


2  =  z  1 1.  a 


s  am  p  1  t  i  o  f  k  o '  z  e  t  a ) 


DOUBLE  N,Nf,N£,H£  !  INTEGERS 

D I M  X  <  S 1 92  > , Y ( 8 1 92 ) , Cos  <  2048 ) 

RED  I M  X  k  0 : Mf -  1 > , Y  <  0 ; Nt - 1 > , Cos  <  0 ; Nf x4 ) 
MRT  X=<0.) 

MFiT  Y=<0.  ) 

T=2. *PlxMf 

FOR  Ns=0  TO  Mf/4 

Cos<Ms)=C0S(:  T*Ms>  !  QUflRTER-CO 

NEXT  Ms 

T  a=  G  p  *  M  u  f  *  M  u  f 

IF  fi=0.  them  220 

Tb  =  Nuf*Muf/  ■vfl*Bs*Bs  ) 

T=FNExp<2. *fl-Ta> 

Tc=T*Ta 
T  d  =  T  a*T  a 

Del f=.  l*Ta/<2. *81  )  !  INCREMENT 

COM  F),  Bs,  2c  ,  Ta,  Tb 

X<0i=. 5»< 1 . -T)  !  TRhPEZOIDR 

FOR  Ns=l  TO  N 

Corr  =  FNKol  <Ns*De1  z)  !  CORRELFITIO 

IF  FiBSCCorrX  1  .  E-38  THEN  340 

Ms  =  Ns  MODULO  Mf  !  COLLFIPSING 

X<M£>=X<Ns)+Corr 

NEXT  Ns 

PRINT  "Final  ualue  of  Corn  =";Corr;" 
MAT  X  =  X*<:Delz) 

CALL  Ffl 14<Mf , Cos<*) , XC«> , Y<*; ) 


QUARTER-COSINE  TABLE 


TRAPEZOIDAL 


!  CORRELATION 
THEN  340 

!  COLLAPSING 


FREQUENCY 


k  o  (  z  e  t  a  >  -  k  1  <  z  e  t  a  > 


Ns  =";Ns 
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3?0  GIMIT 

380  PLOTTER  IS  "GRAPHICS" 

390  GRAPHICS  OH 

400  WINDOW  -4,2,-70,30 

410  LINE  TYPE  3 

420  GRID  1,10 

430  LINE  TYPE  1 

440  FOR  Ni=l  TO  2000 

450  F=Mf*Delf  !  FREQUENCY 

460  W=2.*PI*F 

470  Wl=Tc ''(Td  +  W*W)  !  SHARP  SPECTRAL  COMPONENT 

480  T=Mf*Del2*F 

490  Ns=INT(T> 

500  Fr=T-N£ 

510  W2  =  Fr*X<:N£+l  >  +  <  1  . -Fr  )*>■'<  Ns  )  !  BROAD  SPECTRAL  COMPONENT 

520  PLOT  LGTCF^  ,  10.  ■trLGT':Wl+W2) 

530  NEXT  Ms 

540  Ns=MRXCNs,n 

550  FOR  Ns=Ns  TO  Mt/2 

560  F  =  M£^"(Mt  *Dt  1  z)  !  FREQUENCY 

570  W=2.*PI*F 

580  Wl=Tc  ■  (  T  d  + 1)  *  W  ) 

590  W2=X<N£) 

600  T=W1+W2 

610  IF  T>0.  THEN  640 

620  PENUP 

630  GOTO  650 

640  PLOT  LOTCF-*,  10.*LGT<T) 

650  NEXT  Ms 

660  PENUP 

670  PAUSE 

680  END 

690  ! 

7l0  DEF  FNExpCXini  nus)  !  EXP<-X)  WITHOUT  UNDERFLOW 

710  IF  Xr.,1  riu£>70S.  3  THEN  RETURN  0. 

720  RETURN  EXP  < -Xirh  nus  ) 

730  FNEND 

740  ! 

7‘^f';  DEF  FNKol  (Ze  ta)  !  CORRELATION  k  o  <  zet  a> -t:  1  C  ze  t  a;' 

760  COM  A, Bs, Zc , Ta, Tb 

770  El=FNExo<Zeta)+Zeta-l. 

780  "=Bs*Zeta 

790  F2  =  FNExpCTy  +  T-l  . 

800  R  •,o  =  MAXCO.  ,  1-Zet  a 'Zc  >  f  TRIANGULAR  RHO 

S10  T  =  Ta'*cl+A'*t^.—  R'-  o)— A*R'ho*FNExp(Tb'^^E2) 

820  RETURN  FNE  ■  T  > -FiIF  xp  <  T  a*  <  Ze  t  a- 1  .  > +  2  .  *  A  ) 

830  FNEND 

340  ! 

850  SUB  Fft  14'  DOUBLE  N.REAL  Cos  (  »  )  ,  X  (  ♦  >  ,  Y  <  *■  )  !  N  <  =2"  1  4  =  1  6384  ; 
86c  '  SEE  APPENDIX  A. 3 


0  S  U  B 
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